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ABSTRACT 


A  primitive  equation  spectral  model  using  spherical  harmonics  Is 
formulated  to  study  interactions  between  the  troposphere  and  strato¬ 
sphere  in  association  with  sudden  stratospheric  warmings.  In  order  to 
follow  vertical  wave  propagation  so  important  to  this  process,  the 
model  consists  of  31  levels  with  5  in  the  troposphere  and  the  other  26 
in  the  stratosphere  and  mesosphere.  Using  sigma  coordinates  for  the 
former  and  log-pressure  coordinates  for  the  latter,  separate  equations 
for  each  system  are  combined  to  form  single  matrix  governing  equations. 
The  gradual  introduction  of  planetry  scale  topography  to  an  initially 
balanced  state  representative  of  observed  mean  winter  conditions  in 
the  Northern  Hemisphere  is  used  to  force  changes  in  the  initial  fields 
during  40  day  time  integrations.  Utilizing  the  same  initial  tropo¬ 
spheric  conditions,  three  cases  were  run.  The  first  case  started  with 
a  weak  polar  night  jet  in  the  stratosphere  and  mesosphere.  The  second 
case  had  a  much  stronger  polar  vortex.  For  the  third  case,  a  lid  was 
placed  on  the  troposhere  and  no  interaction  was  allowed  with  the  atmo¬ 
sphere  at  higher  levels. 

Results  of  these  Integrations  indicate  that  realistic  strato¬ 
spheric  warmings  can  be  simulated  by  simple  orographic  forcing*  that 
the  strength  of  the  polar  night  jet  can  be  a  determining  factor  In 
whether  a  warming  becomes  "major"  or  "minor";  that  the  stratosphere 
plays  an  Important  role  In  tropospheric  wave  evolution;  that  the 
indirect  cell  mechanism  is  fundamental  to  the  warming  process;  that 


large  terms  in  the  momentum  and  heat  budgets  tend  to  cancel,  but 
significant  changes  in  the  zonal  wind  and  temperature  do  occur  if 
there  is  a  slight  Imbalance  between  these  terms;  and  that  wave  to 
wave  Interactions  are  an  essential  part  of  sudden  stratospheric 
warmings. 
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CHAPTER  1 


INTRODUCTION 

Since  the  discovery  of  the  sudden  stratospheric  warming  by 
Scherhag  (1952),  this  remarkable  phenomenon  has  been  the  focal  point 
of  numerous  observational  studies  and  more  recently  the  objective  of 
several  modeling  efforts.  Although  our  understanding  of  these  wanning 
occurrences  has  been  greatly  improved  by  these  research  efforts,  many 
questions  still  remain  unanswered.  This  is  particularly  true  of  the 
nature  of  Interactions  between  the  troposphere  and  stratosphere  in 
connection  with  these  warming  regimes.  The  purpose  of  this  thesis  is 
to  explore  this  area  through  numerical  simulation  with  a  sophisticated 
prognostic  model  particularly  well  suited  for  studying  these  processes. 

Sudden  stratospheric  warnings  are  characterized  by  rapid  polar 
temperature  Increases,  sometimes  exceeding  50°K  (within  a  deep  strato¬ 
spheric  level  over  the  period  of  one  week  (Mclntruff,  1978).  Warmings 
have  been  classified  by  the  World  Meteorological  Organization  (WMO)  In 
two  categories.  The  warming  classification  is  major  If  the  temperature 
gradient  and  mean  zonal  wind  direction  reverses  poleward  of  60° 
latitude.  A  minor  warming  Is  one  characterized  by  a  large  temperature 
Increase  anywhere  In  the  stratosphere  with  the  change  occurring  over 
a  short  period,  but  not  meeting  the  conditions  of  a  major  warming.  A 
major  warming  Is  a  sporadic  phenomenon  that  occurs  about  once  every 
other  year.  Minor  warmings  occur  more  frequently  with  sometimes  as 


many  as  three  during  a  winter  season  (Schoeberl,  1978). 

In  absence  of  warming  activity,  the  winter  stratosphere  has  a 
relatively  large  meridional  temperature  gradient  with  cold  polar 
temperatures  due  to  decreasing  ozone  absorption  of  solar  radiation 
during  the  long  polar  night  (Schoeberl,  1978).  As  a  result,  strong 
westerlies  develop  in  response  to  the  thermal  wind  relationship  and 
form  a  cyclonic  vortex  around  the  pole.  The  core  of  these  westerlies 
is  known  as  the  polar  night  jet.  The  jet  maximum  lies  near  the 
stratopause  and  is  actually  located  in  middle  latitudes  further  south 
than  the  boundary  of  the  polar  night  (Craig,  1965).  At  lower  levels, 
the  jet  extends  northward  down  to  the  tropopause.  During  sudden 
warmings,  the  polar  night  jet  undergoes  tremendous  upheaval  and,  as 
in  the  case  of  major  warming,  is  replaced  by  easterlies  over  much  of 
the  polar  stratosphere.  In  some  instances,  as  the  major  warming  of 
1976-1977  (Taylor  and  Perry,  1977;  Quiroz,  1977),  the  warming  and  wind 
reversal  may  extend  all  the  way  to  the  surface  through  the  entire 
troposphere. 

Teweles  and  Finger  (1958)  first  began  to  recognize  from  the  1957- 
1958  major  warming  that  planetary  scale  waves  played  an  important  role 
in  the  warming  process.  This  has  been  verified  by  nearly  every  obser¬ 
vational  study  since  then.  Charney  and  Drazin  (1961)  presented  the 
theoretical  evidence  that  the  only  vertical  propogating  waves  that  can 
penetrate  the  stratosphere  from  the  troposphere  are  those  whose  phase 
speeds  are  westward  relative  to  the  mean  zonal  wind,  but  less  than  a 
critical  velocity  inversely  proportional  to  wavelength.  The  planetary 
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waves,  those  with  zonal  wavenumbers  1-4,  are  usually  the  only  waves 
that  meet  this  criterion  and  then  only  during  the  winter  months  when 
the  upper  atmospheric  flow  is  westerly.  Stratospheric  easterly  flow, 
as  in  the  summer  months  and  as  may  develop  after  a  strong  warming, 
tends  to  trap  most  wave  activity  in  the  troposphere.  According  to 
Houghton  (1978),  if  this  upward  propagation  of  wave  energy  is 
prevented,  radiative  processes  prevail.  After  a  major  warming,  this 
means  a  re-intensification  of  the  polar  night  jet. 

Of  the  planetary  waves,  wave  1  has  the  easiest  time  of  penetrating 
the  winter  stratosphere,  followed  by  wave  2.  Waves  3  and  4  have  a 
much  more  difficult  time  propagating  upward  Into  the  stratosphere,  but 
according  to  Perry  (1967),  wave  3  does  grow  in  the  stratosphere  from 
non-linear  interactions  probably  involving  waves  1  and  2.  Strato¬ 
spheric  warmings  and  their  associated  planetary  wave  activity  have 
been  attributed  to  radiation  effects,  instability,  and/or  forcing  from 
below  (Trenberth,  1973).  However,  all  but  the  latter  appear  to  be 
Insufficient  to  cause  the  large  changes  in  such  a  short  time  scale  as 
with  the  observed  warming.  McIntyre  (1972)  essentially  disproved  the 
instability  as  the  cause  for  the  stratospheric  warmings.  Observed 
changes  in  radiation  during  major  warming  occurrences  (Trenberth, 

1973)  also  tend  to  discount  their  role  in  the  process.  Tropospheric 
forcing,  however,  has  received  support  from  energetics  studies  (Julian 
and  Labitzke,  1965;  Perry,  1967;  Miller  and  Johnson,  1970;  Iwashima, 

1974) ,  observational  studies  (Labitzke,  1977  and  1978;  O'Neill  and 
Taylor,  1979),  and  from  numerical  simulations  (Matsuno,  1971;  Holton, 
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1976;  Schoeberl  and  Strobel ,  1980a  and  1980b). 

Ellassen  and  Palm  (1961)  have  Indicated  that  flow  over  mountains 
produces  gravity  waves,  Inertia-gravity  waves,  and  planetary  waves. 

The  latter  results  from  flow  over  ranges  continental  in  scale  (e.g. 
the  Rockies).  Land  masses  also  radiatively  can  force  large  scale 
wave  patterns.  Because  of  the  greater  topography  and  nonzonal  heat 
contrasts,  the  planetary  waves  1  and  2  are  the  most  important  at  mid¬ 
latitudes  in  the  Northern  Hemisphere  (van  Loon  et  al.,  1973).  On  the 
other  hand,  quasi -stationary  wave  1  and  2  components  cannot  be  readily 
coupled  with  orography  in  mid-latitudes  in  the  Southern  Hemisphere, 
where  the  surface  is  essentially  oceanic  (van  Loon  and  Jenne,  1972). 

It  is  interesting  to  note  that  a  major  sudden  warming  according  to  WMO 
standards  has  never  been  observed  in  the  Southern  Hemisphere  (Schoberl, 
1978),  which  certainly  implies  the  importance  of  topography  in  the 
warming  process,  but  then  begs  the  question  why  are  the  phenomena  so 
sporadic  in  the  Northern  Hemisphere.  Another  unanswered  question  is 
why  pre-1964  major  warmings  were  considered  wave  2  occurrences  and 
every  warming  since  then  has  been  predominantly  associated  with  wave  1 
(Mclntruff,  1978).  Perhaps  this  latter  question  can  be  explained  away 
by  saying  that  the  data  acquisition  by  satellites  and  better  rocket 
and  balloon  observations  have  given  us  a  more  accurate  description  of 
the  phenomena  (Labitzke,  1977),  but  the  doubt  will  remain  unless  in 
the  future  a  strong  wave  2  warming  does  occur.  Atmospheric  models 
(e.g.  Lordi,  1978)  have  easily  produced  wave  2  warmings  adding  to 
this  doubt.  On  the  other  hand,  Labitzke  (1977,  1978)  has  tied  major 
warmings  directly  to  wave  1  enhancement  with  minimal  wave  2  activity 


and  minor  warmings  with  wave  2  enhancement  in  the  troposphere  and 
stratosphere. 

The  sporadic  nature  of  the  warming  is  also  interpreted  in 
several  ways.  One  theory  (Iwashima,  1974)  is  that  the  standing 
components  of  waves  1  and  2  in  themselves  are  not  sufficient  to 
produce  a  warming,  but  if  traveling  components  come  into  phase  with 
the  quasi-stationary  ones  an  explosive  amplification  of  the  wave 
occurs  which  propagates  rapidly  into  the  stratosphere.  Labitzke's 
(1977)  proposal  as  mentioned  above  could  be  another  reason  for  the 
differences  from  season  to  season  and  warming  to  warming.  Perhaps 
wave  2  gets  samplified  more  frequently.  Another  possibility  was 
proposed  by  Koermer  and  Kao  (1980),  where  the  intensity  of  the  warming 
was  inversely  proportional  to  the  strength  of  the  already  established 
polar  night  jet,  which,  if  not  too  strong,  could  reflect  wave  energy. 
Finally,  another  theory  comes  from  Schoeberl  and  Strobel  (1980b), 
who  tie  in  the  orographic  forcing  with  a  stong  tropospheric  jet  as  a 
condition  to  produce  a  major  warming.  All  of  the  above  conditions 
are  feasible  and  could  occur  simultaneously,  since  they  are  not 
contradictory.  This  shows  the  need  for  more  research  into  the 
tropospheric/stratospheric  interaction  field  In  relation  to  strato¬ 
spheric  warmings. 

It  is  generally  accepted  that  as  planetary  waves  propagate  into 
the  stratosphere,  they  set  up  an  Indirect  circulation  in  the  polar 
stratosphere  (Matsuno,  1971).  As  described  by  O'Neill  and  Taylor 
(1979),  this  circulation  consists  of  air  rising  near  the  pole,  moving 
southward  near  the  stratopause,  and  sinking  in  mid-latitudes  with  an 


associated  poleward  heat  flux  at  the  bottom  of  the  cell  in  the  lower 
stratosphere.  Meanwhile,  a  poleward  eddy  momentum  flux  develops  to 
counterbalance  the  coriolis  torque  acting  on  the  southward  movement 
of  air.  Matsuno  (1971)  explained  that  if  there  was  a  sudden  ampli¬ 
fication  of  wave  1  or  2,  this  delicate  balance  would  be  destroyed 
with  net  heating  taking  place  at  the  poles  because  of  increased  heat 
flux  and  the  polar  night  jet  would  weaken  due  to  increased  coriolis 
torque.  Koermer  and  Kao  (1980)  compared  a  major  and  minor  warming 
and  showed  significant  differences  in  their  heat  and  momentum  budgets, 
highlighting  the  importance  of  this  indirect  cell  in  the  warming 
process. 

A  sudden  warming  is  characterized  by  a  downward  propagating  warm 
layer  from  around  45  km  into  the  lower  stratosphere  (Quiroz,  1969, 
1971)  and  tends  to  be  highly  nonzonal  in  nature  (Hirota,  1968).  The 
magnitude  of  the  effects  of  the  warming  process  are  emphasized  by  not 
only  a  nearly  compensating  cooling  in  the  mesosphere  (Labitzke,  1972; 
Hirota  and  Barnett,  1977),  but  also  by  a  weaker  cooling  in  the  lower 
equatorial  stratosphere  of  about  2°  for  every  30°  to  40°  of  warming 
that  takes  place  in  polar  latitudes  (Houghton,  1978).  All  of  these 
effects  make  a  highly  dynamic  system  out  of  the  stratosphere,  which 
was  once  believed  to  be  a  docile  part  of  the  atmosphere. 

The  best  known  simulations  of  stratospheric  warmings  were  made 
by  Matsuno  (1971)  and  Holton  (1976).  Matsuno  used  a  quasi -geostrophic 
model.  It  was  a  mechanistic  model  with  forcing  specified  at  the 
bottom  boundary  of  the  model  which  corresponded  to  the  10  km  level . 
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The  troposphere  was  not  included  in  this  formulation.  The  forcing 
consisted  of  a  single  zonal  harmonic  component.  Interaction  was 
allowed  only  between  this  wavenumber  and  the  mean  flow.  His  model 
produced  many  of  the  basic  characteristics  of  observed  warmings. 

Quiroz  et  al .  (1975)  compared  this  model  to  recently  available 
satellite  data  that  generally  supported  the  model  results.  They 
did  note  several  discrepancies,  such  as  cooling  did  not  occur  in 
the  upper  stratosphere  in  Matsuno's  model  as  has  been  observed; 
the  simulated  westerlies  In  the  mesosphere  were  much  stronger  than 
observed;  and  the  changes  during  the  first  1-10  days  are  too 
monotonic  in  comparison  to  observed  vacillations.  According  to 
Schoeberl  (1978),  part  of  Matsuno's  problem  was  due  to  the  fact 
that  neither  damping  nor  friction  was  Included  in  his  calculations. 

This  did  not  allow  the  mean  zonal  flow  to  relax  after  eddy  forcing 
and  also  caused  the  critical  level  to  numerically  reflect  any  wave 
advance. 

Holton  (1976)  used  a  primitive  equation  semi-spectral  model. 

This  model  was  also  linear  and  confined  to  the  upper  atmosphere. 

The  forcing,  like  Matsuno's,  was  turned  on  slowly  at  the  tropopause 
and  was  held  steady  after  reaching  a  maximum  after  several  days  of 
integration.  His  results  generally  tended  to  confirm  many  aspects  of 
Matsuno's  model.  Holton  did  use  damping  and,  perhaps  as  a  consequence, 
his  model  did  not  show  the  downward  propogation  of  the  zero  wind  line 
nor  did  he  get  the  double  warming  In  the  stratosphere.  Quiroz  (1969) 
has  basically  discounted  the  presence  of  downward  propagating  critical 
levels,  but  Quiroz  et  al.  (1975)  believe  a  double  warming  as  simulated 


by  Matsuno  is  possible.  Holton's  easterlies  are  much  more  conserva¬ 
tive  than  Matsuno's;  but  in  the  case  of  Holton's  wave  1  forcing,  the 
wind  reversal  did  not  extend  deep  enough  into  the  lower  stratosphere 
to  be  classified  a  major  warming  by  WMO  standards. 

Using  a  high  vertical  resolution  spectral  model,  Lordi  (1978) 
went  further  and  simulated  sudden  warmings  allowing  not  only  for  wave 
and  mean  flow  but  also  for  wave-wave  interactions  deemed  so  important 
by  Perry  (1967).  Starting  with  essentially  the  same  initial  condi¬ 
tions  as  Matsuno  (1971)  and  Holton  (1976)  and  employing  the  same 
single  zonal  harmonic  forcing,  he  was  able  to  produce  sudden  warmings. 
Using  a  linearized  version  of  his  model,  Lordi  was  also  able  to 
simulate  warmings,  but  the  contrast  between  these  results  and  those 
obtained  using  the  "non-linear"  version  was  quite  remarkable.  In  this 
thesis,  the  term  "linear"  means  Interaction  is  allowed  only  between 
the  wave  and  mean  flow.  "Non-linear"  means  the  wave  can  interact  with 
other  waves  and  the  mean  flow.  For  wave  1  forcing,  the  non-linear 
simulation  resulted  in  a  more  rapid  warming  than  the  linear  case  and 
Indicated  the  contribution  of  higher  harmonics  to  the  warming  process. 
For  wave  2  forcing,  the  opposite  was  true  as  the  linear  case  was  more 
dynamic.  Perhaps  this  was  due  to  the  fact  that  since  only  wave  2  was 
forced,  wave  1  remains  unaffected  and  cannot  play  an  essential  role 
in  the  warming  process.  Hence,  wave  energy  was  only  transferred  back 
In  the  mean  flow  or  to  wave  4  (the  point  of  truncation),  which  may 
detract  from  wave  2 ' s  ability  to  produce  a  strong  warming.  In  the 
linear  case,  wave  2  does  not  yield  any  energy  to  higher  harmonics  and 
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may  cause  the  simulated  warming  to  progress  more  strongly.  These 
results  clearly  Indicate  that  both  waves  1  and  2  must  play  an  impor¬ 
tant  role  in  the  warming  process.  The  evolution  of  Lordi's  wave  1 
warming,  in  particular,  which  reflects  wave  2  contributions  as  well, 
can  be  compared  closely  with  the  evolution  of  the  warming  of  January 
1973  as  depicted  by  Mclnturff  (1978).  Another  key  observational 
consistency,  which  is  not  clearly  evident  in  the  data  presented  by 
other  researchers  mentioned  above,  is  the  evolution  of  the  secondary 
jet  maximum  in  the  westerlies  around  35  km  at  higher  latitudes  during 
the  pre-warming  stage  as  noted  by  Quiroz  et  al .  (1975).  This  is  pro¬ 
duced  in  all  of  Lordi's  Integrations  and  precedes  the  sudden  warming. 
The  author  has  made  detailed  observational  studies  of  the  lower 
stratosphere  during  the  major  warmings  of  1973-1974  and  1976-1977. 
These  essentially  wave  1  warmings  are  in  close  agreement  with  Lordi's 
wave  1  non-linear  results  in  that  there  is  no  evidence  of  strong 
easterly  development  In  the  mid-latitudes  of  the  lower  stratosphere  as 
evolved  in  wave  2  forcing  cases  and  to  some  extent  the  wave  1  linear 
case.  In  all  of  Lordi's  runs,  the  stratospheric  warmings  were  accom¬ 
panied  by  cooling  in  the  upper  stratosphere  and  mesosphere  as  well  as 
cooling  in  the  equatorial  stratosphere. 

The  most  recent  simulated  stratospheric  warming  calculations  were 
done  by  Schoeberl  and  Stroberl  (1980a).  They  used  a  model  similar  to 
that  of  Matsuno  (1971),  but  Included  solar  radiative  heating  due  to 
O2  and  O3  absorption,  Newtonian  cooling  and  Rayleigh  friction.  They 
again  forced  with  a  single  zonal  harmonic  and  compared  evolution  of 
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their  calculations  with  and  without  damping.  Without  damping,  they 
obtained  results  similar  to  Matsuno's  (1971)  with  strong  easterlies 
and  fairly  rapid  reversals.  With  large  damping,  the  warmings,  as 
should  be  expected,  occurred  at  a  much  slower  rate  and  the  easterly 
maximum  that  did  develop  was  much  less  than  those  obtained  in  the 
undamped  case.  They  indicated  that  damping  controls  the  amplitude 
that  a  warming  can  achieve.  They  also  reinforced  the  concept  of 
Holton  (1976)  that  lateral  critical  levels  do  play  an  essential  role 
in  the  warming  process.  Their  wave  2  warming  evolved  more  rapidly 
than  their  wave  1  warming  in  contrast  to  the  results  of  Lordi  (1978), 
who  used  a  fully  non-linear  model.  Schoeberl  and  Strobel  (1980a) 
also  mply  that  warmings  require  strong  westerlies  in  the  lower 
stratosphere  in  contrast  to  the  observational  data  presented  by 
Koermer  and  Kao  (1980). 

Schoeberl  and  Strobel  (1980b)  recently  added  a  troposphere  to 
their  previous  model  and  orographically  forced  waves  in  the  troposphere. 
Although  in  some  of  their  cases,  waves  1  and  2  were  forced  simultan¬ 
eously,  the  waves  were  allowed  to  interact  with  each  other  through  the 
mean  flow  only.  Topography  was  gradually  introduced  over  a  15  day 
period.  Because  of  the  large  damping  used,  it  took  85  days  to  develop 
a  major  warming  and  then  only  after  the  tropospheric  jet  was  increased 
to  48  ms-1  from  36  ms-*  which  they  used  in  a  prevous  attempt  and 
failed  to  produce  a  warming.  Although  they  presented  only  a  brief 
account  of  their  results,  they  concluded  that  wave  1  is  much  slower 
developing  than  wave  2  and  therefore  major  warmings  must  be  wave  2 


phenomena,  contradicting  the  observational  data  of  Labltzke  (1977, 
1978).  They  also  Imply  that  a  strong  tropospheric  jet  and  a  strong 
polar  night  jet  are  essential  for  a  strong  warming,  conditions  which 
are  not  supported  by  Koermer  and  Kao  (1980).  Perhaps  the  reasons  for 
these  apparent  discrepancies  are  that  the  damping  used  by  Schoeberl 
and  Stroberl  (1980b)  was  too  large  and  secondly  wave-wave  Interactions 
were  not  permitted.  It  certainly  does  not  appear  to  be  realistic  that 
sudden  warmings  should  take  so  long  to  develop. 

Perry  (1967)  stressed  not  only  the  Importance  of  non-linear 
interactions  In  the  stratosphere,  but  also  the  importance  of  wave-wave 
interactions  In  the  troposphere  as  part  of  the  sudden  warming  process. 
So  far  In  a  detailed  model  designed  to  study  the  warming  phenomena 
only  Lordl  (1978)  has  allowed  for  wave-wave  interactions  In  the 
stratosphere,  but  the  tropospheric  contribution  was  limited  to 
continuous  forcing  applied  at  the  tropopause.  This  forcing  poses  an 
artificial  condition  on  the  system  as  evident  from  the  latitude  time 
sections  of  geopotential  for  waves  1  and  2  presented  by  O'Neill  and 
Taylor  (1979).  The  model  of  Schoeberl  and  Strobel  (1980b)  Introduced 
a  troposphere  to  replace  the  continuous  forcing  at  the  tropopause  with 
a  more  realistic  generation  of  wave  1  and  2  by  orography.  However, 
this  model  does  not  allow  for  wave-wave  interactions,  which  from 
observational  data  Perry  (1967)  had  found  so  Important.  These 
conditions  may  explain  why  none  of  the  above  models  have  been  able  to 
produce  warmings  with  wind  reversals  comparable  In  vertical  depth  as 
that  reported  for  the  winter  of  1976-1977  (Taylor  and  Perry,  1977). 


The  purpose  of  this  thesis  then  is  to  simulate  a  sudden  strato¬ 
spheric  warming  with  a  completely  non-linear  model  which  also  includes 
a  troposphere.  The  model  will  be  forced  by  a  combination  of  zonal 
harmonic  components  1  and  2  by  using  topography  to  simplistically 
represent  the  gross  orographic  features  of  the  Northern  Hemisphere. 

This  allows  for  variable  forcing  of  the  stratosphere  from  the  tropo¬ 
sphere.  Unlike  many  general  circulations  models,  this  model  will 
include  sufficient  vertical  resolution  to  closely  follow  vertical  wave 
propagation  that  Is  so  Important  in  the  warming  process  and  in  any 
tropospheric-stratospheric  Interactions. 

In  this  undertaking,  we  will  use  a  modified  version  of  Lordi's 
(1978)  model  for  the  stratosphere  and  mesosphere.  It  will  be  coupled 
with  a  tropospheric  spherical  harmonic  model  adapted  from  the  work  of 
Bourke  (1974)  and  Hoskins  and  Simmons  (1975)  with  modifications. 

Studying  only  planetary  scale  waves,  we  then  want  to  examine  the  role 
of  orography  in  a  warming  event;  the  effect  of  the  strength  of  the 
polar  night  jet;  the  pre-eminence  of  wave  1  or  2  in  the  warming  process; 
importance  of  an  upper  atmospheric  representation  in  long  range 
numerical  simulations;  and  the  feedback  effects  that  the  stratosphere 
may  have  on  the  troposphere.  These  latter  effects  may  contribute  to 
long-term  tropospheric  climatic  fluctuations  as  suggested  by  Kasahara 
et  al .  (1973). 


CHAPTER  2 


MODEL  FORMULATION 

The  model  to  be  formulated  in  this  chapter  basically  consists  of 
two  parts.  The  upper  part  of  the  model  includes  the  stratosphere  and 
mesosphere  and  the  basic  equation  development  follows  very  closely 
after  Lordi  (1978).  The  lower  tropospheric  part  of  the  model  is 
loosely  patterned  after  the  tropospheric  models  of  Bourke  (1974)  and 
Hoskins  and  Simmons  (1975).  The  governing  equations  for  each  section 
of  the  model,  although  similar  in  many  respects,  are  also  radically 
different  primarily  because  of  the  different  vertical  coordinate 
systems  employed  in  each  part.  However,  these  somewhat  independent 
equations  can  be  merged  into  common  matrix  equations  which  allow  us 
to  perform  semi-implicit  time  differencing  on  the  entire  system  at 
the  same  time. 


2.1  Vertical  Coordinate  Systems 

In  his  model,  Lordi  used  a  log-pressure  coordinate  system  defined 
as  follows: 

z  =  -  H  in  (p/ps>,  (1) 

where  z  is  the  height  coordinate;  p  is  the  pressure;  p^  is  the  mean 
surface  pressure  taken  to  be  1016  mb,  and  H  is  the  scale  height 
defined  by 
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H  =  RTm/g  ,  (2) 

where  R  is  specific  gas  constant  of  dry  air,  g  the  gravitaional  acceler¬ 
ation  g,  and  Tm  a  mean  atmospheric  temperature  assumed  to  be  244°K. 
While  this  system  is  satisfactory  for  the  upper  atmosphere,  a  modified 
sigma  coordinate  system,  first  introduced  by  Phillips  (1957),  is  used 
for  the  troposphere  in  order  to  handle  orography  which  is  so  important 
to  any  tropospheric  model. 

This  0-system  is  defined  as  follows: 

a  =(p  -  PR)/(P*  -  PR)  ,  (3) 

where  pr  is  a  reference  pressure  indicating  the  upper  boundary  of  the 
tropospheric  model  and  lower  boundary  of  the  upper  model;  p^  is  the 
surface  pressure;  and  p  is  the  pressure  at  a  given  level.  It  follows 
from  (3)  that  a  =  0  at  p  =  p^  and  a  =  1  at  p  =  p^.  In  order  to  deter¬ 
mine  the  reference  pressure  level,  we  can  solve  (1)  for  p  to  obtain 

p  =  ps  exp  (-  z/H)  .  (4) 

In  his  model,  Lordi  had  a  lower  boundary  height  of  zb  =  12  km.  Hence, 
from  (4)  we  can  compute  a  reference  pressure  pr  *  189.5  mb,  which  will 
make  the  o  =  0  level  of  the  troposphere  coincide  exactly  with  the  12 
km  level  of  the  upper  atmospheric  model. 

Figure  1  depicts  the  vertical  grid  structure  of  the  coupled 
models.  In  the  stratosphere  and  mesosphere,  az  =  3  km  with  26 
prognostic  levels,  which  will  generally  be  represented  with  the  index 
K.  The  index  r  is  used  to  indicate  the  5  prognostic  levels  of  the 
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troposphere  where  a o  =  .2.  Vertical  velocity  (w  or  o)  and  geopotential 
W  are  diagnostic  variables  and  values  are  computed  for  diagnostic 
levels  centered  between  the  prognostic  levels.  Vorticity  (c), 
divergence  (0)  and  temperature  (T)  are  the  prognostic  variables  common 
to  both  models.  Additionally,  log  pressure  (q)  is  a  prognostic  vari¬ 
able  in  the  troposphere,  for  boundary  conditions,  we  assume  vertical 
velocity  is  zero  at  the  top  of  the  upper  model  and  the  bottom  of  the 
lower  model.  We  also  specify  surface  geopotential  as  a  bottom 
boundary  condition.  Vertical  derivatives  at  the  bottom  of  the  tropo¬ 
sphere  will  be  specified  in  the  parameterization  of  vertical  diffusion 
to  be  discussed  later  in  this  chapter. 

2.2  Governing  Equations 

The  momentum,  thermodynamics,  continuity,  and  hydrostatic 
equations  equations  serve  as  the  foundation  for  both  models.  Although 
the  development  of  the  upper  atmospheric  model  equations  will  not  be 
discussed  in  detail,  mention  will  be  made  of  slight  modifications  and 
essential  equations  as  developed  by  Lordi  (1978)  will  be  introduced  to 
complete  the  formulation.  In  the  remaining  part  of  this  chapter,  we 
will  develop  the  tropospheric  governing  equations,  give  the  strato¬ 
spheric  counterpart,  form  common  matirx  equations,  and  introduce  the 
time  differencing  schemes  use  to  perform  model  integrations. 

After  Kasahara  (1974),  the  basic  equations  for  a  dry  atmosphere 
in  o-coordi nates  can  be  written  as  follows: 

+ 

dv  =  -  fk  X  V  -  v*  -  qgl  VP  +  F  ; 

TTf  oP~+Pr  * 


(5) 
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dT  =  RJ  p  +  Q_  ;  (6) 

"eft  cp  ( dP*  +pr  )  cp 


3P  > 

_1_  +  V  •  (P  V)  +  P  l£  =  0  (7) 

8t  *  *3 o 


RIP 

34  =  -  * 

So"  OP^+PR  ’ 


(8) 


where  V  is  the  horizontal  verlocity;  f  is  the  coriolis  parameter;  k 

* 

is  the  vertical  unit  vector;  v  is  the  horizontal  gradient  operator; 

♦ 

*  Is  the  geopotential.  T  Is  the  temperature;  P^  =  p#  -  pp;  F  is  the 
horizontal  frictional  force;  cp  is  the  specific  heat  of  dry  air  at 
constant  pressure;  p  h  dp/dt,  the  vertical  velocity  in  pressure 
coordinates;  Q  is  the  diabatic  heating;  o  i  do/dt  is  the  vertical 
velocity  in  a-coordi nates;  and  t  represents  time. 

Using  equation  (5),  we  will  be  able  to  derive  both  vorticity  and 
divergence  tendency  equations.  Those  derivations  along  with  manipu¬ 
lated  versions  of  ( 6 ) - ( 8 }  will  serve  as  the  basic  model  equations.  In 
the  following  sections,  we  will  develop  these  equations. 

2.2.1  Vorticity  Equations 

Expanding  the  total  time  derivative  of  (5)  and  rearranging  terms, 
the  momentum  equation  can  be  written  in  the  following  form: 


3V/3t  =  -  (c+f)k  Xi-V(M  V*V)  -  RTVq  -  o(3V/3o)  +  F 


(9) 


¥ 

l 


where  f  is  the  coriolis  parameter,  q  =  *n(oP  +p  )  and  ;  i  k  •  V  x  V  , 

*  R 

the  vertical  component  of  relative  vorticity.  If  we  let  T  =  T(o)  + 
T',  (i.e.,  temperature  is  composed  of  a  layer  mean  and  a  deviation) 
and  take  k  •  V  x  (9),  it  follows  that 


3;/at  =  -  V  •  (c+f )V  -  k  .  V  x  ( RT 1 Vq+  o  3V/3o  -F)  .  (10) 

•+ 

By  the  Hemholtz  theorem  V  can  be  expressed  as  the  sum  of  non- 
divergent  and  irrotational  parts  as  follows: 

(in 

where  ip  is  the  stream  function  and  *  is  the  velocity  potential.  Hence, 
it  follows  that 

c  =  k  •  v  x  V  =  V2*  ,  (12) 

D  =  V  >  V  =  V2X  ,  (13) 

where  D  is  the  horizontal  divergence.  In  spherical  coordinates,  we 
can  express  (11)  in  component  form 

U  =  -  J_  +  1  3X  , 

a  7$  a  cos  $  T\  (14) 

v  =  1  ^  +  1  3* 

a  cos  $  3X  a  3<f>  (15) 


where  u  Is  the  zonal  velocity  component;  v  is  the  meridional  velocity 
component,  a  is  the  earth's  radius;  $  is  the  latitude;  and  x  is  the 
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longitude.  If  we  further  define 


U  =  u  cos  tp  , 

(16) 

V  i  v  cos  , 

(17) 

it  follows  from  (14)  and  (15)  that 

u  =  -  COS  4)  +  l  3x. 

a  3  $  a  3  x  (18) 

V  =  1. +  cos  <p  3  x 

a  3x  a  3  <j>  (19) 

-► 

We  can  also  express  the  frictional  force  F  in  component  form  as 

F  »  FAi  +  Fj  (20) 

where  F\  and  F ^  are  the  longitudinal  and  meridional  components, 
respectively.  The  parameterization  of  these  components  will  be 
detailed  later  in  this  chapter. 

In  order  to  control  spectral  blocking  (Puri  and  Bourke,  1974) 
that  results  from  horizontal  truncation,  diffusion  is  applied  to  the 
vorticity  equation.  Similar  to  Lordi  (1978),  we  used  the  linear 
diffusion  mechanism  of  Leith  (1971).  This  technique  damps  deviations 
from  equilibrium  fields  and  not  the  equilibrium  fields.  With  this 
specification  added  and  using  (12),  ( 1 4 ) - ( 20 ) ,  we  can  express  the 
vorticity  equation  (10)  in  spherical  coordinate  form  as  follows: 

*  -  1  .  [j£  +  COS  »  ID  +  K  iV2(;-c  )  + 

3t  3 x  3$  h  e  2 

a  cos  4  a 


(21) 
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where  we  define 

A  h  U  +  f)U  +  o ( 3V/3o )  +  (RT'/a)cos  4»(3q/S$)  -  F$  cos  *  (22) 

B  h  U  +  f)V  -  o(3U/3o)  -  (RT7a)(3q/3X)  +  Fx  cos  *  .  (23) 

is  the  horizontal  diffusion  coefficient.  The  e  subscript  represents 
equilibrium  vorticity  at  the  start  of  model  integrations. 

The  parallel  stratospheric  equation  as  given  by  Lordi  (1978)  is 
as  follows: 

2 1  j 

W2!’  =  -  1  fSfi  +  cos  <I>  +  K  |v2(c-c  )  +  ]  (24) 

3t  a  "cos?  <j>  3X  3<f  h  e  ^2 

where 

=  U  +  f)U  +  w( 3V/3z)  -  F ^  cos^>  ,  (25) 

t  =  U  +  f)V  -  w(3'J/3z)  +  Fx  cos^  ,  (26) 

for  w  =  dz/dt  and  z  is  defined  by  (1). 

Combining  (21)  and  (24),  into  a  single  expression,  we  can  write 

3c/3t  =  Z  +  K h[V2(c-Ce)  +  2(C  -  Ce)/a2]  (27) 

where  Z  take  on  the  values  of  the  non-diffusive  terms  on  the  right 
side  of  (21)  for  the  troposphere  and  the  right  hand  side  of  (24)  for 
the  stratosphere  and  mesosphere. 


2.2.2  Divergence  Equation 

The  divergence  tendency  equation  can  be  obtained  by  taking 
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(9) ,  which  yields 


S(V-V)  =  V  *  [-  (c  +  f)k  x  V]  -  ?2(4»  +  V-V  +  RT) 
51  T"  q 


-  V  •[RT'Vq  +  a  (3V/3a)  -  F] 


(28) 


Expanding  (28)  in  spherical  coordinates,  introducing  horizontal 
diffusion,  and  grouping  terms  as  with  the  vorticity  equation,  we 
obtain 
2 

I!.*  =  ..V  .—  [  jjj.  -  cos  4>  M.]  -  v2(E^+RTq)+K.  [72D+2(D/a2)]  (29) 

a  cos^  dA  3<f  n 

where 

♦  ♦  2  2 

E  =  V-V  =  U  +  V  .  (30) 

"T-  2  col?"  * 

For  the  upper  portion  of  the  model,  we  have 

2  2  2 

9V  x  =  1  [i|  -  cos  «  jAj  -  v  (E+*)  +  K  [V  D  +  Z(D  )].  (31) 

3t  a  cos^  $  ax  a<{>  ^ 

As  a  single  combined  equation,  we  can  write  the  divergence  tendency 
equation  as 

3D/3t  =  V  -  +  cRTq)  +  Kh[V2D  +  2( D/a2 ) ]  ,  (32) 


where  c  =  0  for  the  upper  portion  of  the  model  and  c  =  1  for  the 
troposphere  and  where  V  takes  on  the  remaining  values  on  the  right  of 
(29)  that  are  not  explicitly  shown  in  (32)  and  on  the  right  of  (31) 
in  a  similar  fashion  for  the  lower  and  upper  portions  of  the  model 


respectively . 


2.2.3.  Frictional  Force  F 

In  this  section  we  will  define  explicitly  the  frictional  cornpon- 
-► 

ents  F  and  F  of  F  that  were  alluded  to  in  the  development  of  the 

x  <f> 

vorticity  and  divergence  tendency  equations.  For  the  upper  portion 

of  the  model  F  and  F  are  characterized  by  Rayleigh  friction.  In 

X  <)> 

the  troposphere,  vertical  diffusion  is  also  parameterized  in  these 
components  along  with  Rayleigh  friction. 

Following  the  vertical  stress  formulations  of  Bourke  et  al . 
(1977)  and  the  Rayleigh  friction  method  used  by  Lordi  (1978),  we  can 
express  the  tropospheric  components  of  F  as  follows: 


3t 

F  =  (9_)  _A  -  F  (u-u  )  , 

X  p  3o  re 


(33) 


3t 

F  =  (9_)  _ i  -  F  v  ,  (34) 

4>  p  3o  r 
★ 


where  F  is  the  Rayleigh  friction  coefficient,  u  is  the  initial 
r  e 

equilibrium  zonal  wind  components;  and 


T 

X 


2 

=  P  (9_)K  , 

p  v  3o 
r* 


(35) 


T 

$ 


2 

=  p  (g_)x  *1 , 

p  V  3o 


(36) 


for  density  p  and  vertical  diffusion  coefficient  K  defined  as 

v 


22 


K  =  p(9_>2|3V|  ,  (37) 

v  |7o| 


where  v  Is  the  mixing  length  which  we  assume  has  a  value  of  30  meters 
for  a  >  .1  and  is  zero  for  o  <  .1. 


The  lower  boundary  specification  is  given  by 


t  =  -  p  C  I V  I  u 

*\  H  d  N  N  (38) 


T  =  -  p  C  I V  I  V  , 

*4>  N  d  N  N  (39) 


where  the  N  subscript  denotes  the  lowest  prognostic  level, 
drag  coefficient  assumed  to  be  .0025. 


C  i s  the 
d 


2.2.4  Vertical  Velocity  Expressions 

Before  we  can  proceed  with  the  development  of  the  remaining 
equations,  we  must  first  derive  diagnostic  expressions  for  vertical 
velocities  of  the  various  coordinate  systems  encountered  and  the 
relationships  between  them.  Starting  with  equation  (3),  it  follows 
that 

o  =  i_  [p  -  o  _ 1  ]  .  (40) 

P*  3t 


From  (1 ),  we  also  have 

w  =  -  ( H/p ) p  . 

Since  at  p  =  p  ,  a  =  0  we  can  combine  (40)  and  (41)  to  obtain 
R 


(41) 
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0  =  -  *_  P 

r  nr  r 

★ 


where  the  R  subscript  represents  reference  pressure  level  evaluation. 

Integration  of  the  continuity  equation  (7)  with  respect  to  o 
from  0  to  a  using  (42)  yields 

3P*  a  w  n 

a  —  +  /  V  •  (P.V)da  +  P  o  +  R  -R-  =  0  .  <43) 

3t  *  *  H 

Extending  the  integration  of  (43)  to  the  earth's  surface  and  using  the 

bottom  boundary  condition  that  a  =0,  then  we  have 

o=l 


3P*  1  w  D 

—  =-/f  .(P>o  .Jgl- 


Substitution  of  (44)  into  (43)  yields 


w  p  o 

r\  . 


o[-J  v  •  (P^V)do  -  R  Rl  +  /  7  •  (P+V)do  +  P^o  +  R  R  =  0  (45) 


and  solving  (45)  for  o,  we  can  write 


*  1+  -r  4  W  D 

o  ■  «_  /  V  •  (P  V)do  -  j  7  •  (P  V)do  +  (o-l)  RR.  , 
P*  o  *  P*  o  * 


which  gives  us  a  diagnostic  expression  for  a. 

If  we  substitute  (46)  into  (40)  and  solve  for  p,  we  obtain 


•  1  +  +  a  *  -*■  W  P  dP 

p  =  o  /  7  •  (P  V)do  -  /  7  •  (P  V)do+(o-l)  ft  ft  +  o  *  . 

0  *  0  *  H^r 


(47) 
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Si  nee 


then  substitution  of  (44)  into  (48)  yields 

dP*  1  ♦  +  w  p  ♦  + 

—  =  -  /  v  .  (P*V)do  -  +  V  •  VP,  • 

dt  o  “ 

Using  (49)  in  conjunction  with  (47)  yields 


.  o  +  +  w  p  +  + 

p  =  -  /  V  •  (P^V)do  -  R  R  +  oV  •  VP*  , 


(48) 


(49) 


(50) 


a  diagnostic  expression  for  p. 

Although  tropospheric  vertical  velocity  computations  in  height 
coordinates  are  not  required  in  any  model  calculations,  it  is  computed 
for  output  and  display  ourposes.  Tropospheric  vertical  velocities  can 
then  be  related  more  meaningfully  with  the  vertical  velocities  of  the 
upper  portion  of  the  model,  even  though  those  values  are  not  with 
respect  to  a  true  height  system.  From  the  relationship  that  w(3p/3z) 

=  o( 3p/3a) ,  we  can  solve  for  w  and  substitute  the  hydrostatic  equation 
(8)  to  obtain 


-RTP,  •  _ 

g(oP*+pR)° 


(51) 


Since  we  need  only  approximate  values  for  output  purposes,  we  assume 


T  =  T  . 
m 


I 


25 

In  the  upper  part  of  the  model,  vertical  velocity  Is  expressed  In 
two  ways  as  follows: 


ZT0P  +  + 

w *  =  /  (exp(-z)  7  •  V jdz  (52) 

Z  " 


w  =  w*  exp(Z)  .  (53) 

H 


Therefore,  at  the  lower  boundary  at  the  Interface,  (52)  and  (53)  imply 

zn  ZT0P  +  > 

w  =  exp (_5.)  J  (exp(-Z)v  •  V jdz  ,  (54) 

R  Hz 

R 

which  can  be  used  in  conjunction  with  (46)  and  (50)  in  isolating  the 
divergence  in  the  thermodynamic  and  continuity  equation  in  preparation 
for  semi-implicit  time  differencing. 


2.2.5  Thermodynamic  Equation 

Expanding  the  thermodynamic  equation  (6)  and  expressing  tempera¬ 
ture  in  terms  of  a  layer  mean  and  deviation  allows  us  to  write 


3T‘  3  _  y  .  vT'  -  o  a(T+T')  +  R(T+T')  J  +  Q__  .  (55) 

3t  3 o  C  (oP  +p  )  C 

p  *  R  p 


Further  expansion  of  the  gradient  operator  v  and  writing  (55)  in 
spherical  coordinates,  it  follows  that 


33‘  =  '  V  +  COS  $  d-VT'  1  +  DT' 

TT  2  L^1T  T* 


a  cos  $ 


o  a(T+T‘)  +  R(T+T‘)  p  +  Q_  . 

3a  C  ( aP  +p  V  C 
p  *  R  p 


Substitutions  of  (46)  for  a  for  the  mean  vertical  gradient  part  of  (56) 
and  (50)  for  p  and  expanding  the  horizontal  gradient  operator  in  these 


expressions,  i.e,  v  •  (P^V )  =  P^D  +  V  •  vp^,  we  can  write 


of  =  _  1  [3UT1  +  cos  4>  SVT'l  +  DT*  +  Q_  -  a  ill 

Tt~  2  3X  34  c  3a 

a  cos  4>  p 


1  O 

(3T/3o)[|  J  (V  •  VPjdo  -  ~  /  (V  •  VPjdo] 
*  0  y*  0 


.  0  ♦  -*•  *  w  pn 

RT  [  /  (P.D  +  V  •  VP  Ida  -  oV  •  VP  +  R  R  ] 

cp(oP/PR)  0  H 


_ AT _  [  /  (V  •  VP  )do  -  aV  •  VP  ] 

C  (aP  +p  )  0  *  * 

p  *  R 


_  1  a  w  p 

H  [a  /  Oda  -  /  Ddo  +  (a  -1)  R  R 
3o  o  0  PH 


RTP  a  w  p 

—  — *  —  [  /  Dda  +  44  ] 
c  (aP  +p  )  A  V  H 

p  *  R  * 


Similar  to  Lordi  (1978),  we  express  the  diabatic  heating  as 
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Q/C  =  a( z ) ( T 1  -  T')  (58) 

P  e 

which  is  based  on  the  work  of  Leovy  (1964),  where  T'  is  the  equili- 

e 

brium  temperature  deviation  and  a(z)  is  the  Newtonian  heating-cooling 
coefficient  given  by  Holton  (1976): 

a(z)  =  [1.5  +  tanh(Z -35000 Q1  x  lo'V*  ,  (59) 

7U00. 

where  the  height  corresponding  to  the  a-coordinate  level  is  computed 
from  mean  and  deviation  geopotential  fields. 

In  order  to  have  divergence  involved  with  mean  quantities,  we 
define 


6  =  6(a)  +  6* 


oP^+Pr 


(60) 


5  =  T(o)  +  6’  =  (1/PJ  .  (61) 

Using  (54),  (60),  and  (61)  to  isolate  divergence  with  mean  terms  and 
adding  horizontal  and  vertical  diffusion,  we  can  express  (57)  as 
follows: 


I 
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3T‘  =  -  1  PUT1  +  cos  *  aVT’l  +  Dr1  +  afzKT'-T') 

7T~  2  3A  3$  e 

a  cos  4> 


+  K  [V  (T'-T‘ )]  +  F  -  a(3T73o) 
h  e  v  T 


-  1++  0  +  +  w  p 

(3T/3o)[£  /  ( V •  VP  )da  -  1_  /  { V *VP  )do  +  6'(o-l)_0] 
P  0  *  P  0  *  H 


a  +  +  +  •+•  w  p 

RT_  [  /  (P  D+V *VP  )da  -  aV*VP  +  R  R  ] 
C  ( aP  +p  )  0  *  *  *  H 


RT 


w  p 


[  /  ( V*VP  }do  -  oV *7P  +  P  6'  R  R1 

c(oP+p)o  *  *  *  H 

P  *  R 

_  o  w  p 

-  RI  3*C  /  Dda  +  T  JL2  ] 


cp  o 


H 


1  a  n  7  T0P 

-  (3T/3o)[a/  Dda  -  /  Dda  +  T (a-1  )lB.  exp(JL)  /  D  exp(zi)dz] 

oo  H  H  z  H 

R 

a  nZ  Z  TOP 

-  RT  0  [  /  Dda  +  T  lR  exp(_R)  /  D  exp(il)dz]  ,  (62) 

co  H  H  z  H 

P  R 

where  F  denotes  the  vertical  diffusion  term.  After  Bourke  et  al . 
v  T 

(1977),  F  is  defined  as  follows: 
v  T 


Ft  =  (g/p  )0n/3a)  , 
v  T  * 


(63) 


where 


n  =  lip  (g/p  )K  ( 3 8/3 a) 
*  v 


l 


(64) 
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Here  K  is  defined  by  (37) ;  0  =  T(p  /p)  P  is  the  potential  tempera- 

R/c 

ture  with  p  =  1000  mb,  and  n  =  (p/p  )  P.  For  the  lower  boundary 
o  o 

specification  we  assume  thermal  insulation,  i.e.  =  0. 

For  the  statosphere,  we  have  a  thermodynamic  equation  of  the 
following  form: 


!£•  =  -  1  [  a(UT')  ♦  cos  ♦  a(VT‘)i  +  dt*  -  w  |H 

dt  2  »a  d<f>  Ji 

a  cos  $ 


2  2  TOP 

+  K  [V  (T'-T1)]  -  iLH  +  a(T‘-T‘)  -  exp(£)  j  D  exp(l£)dz,  (65) 
h  e  CpH  e  R  H  z  H 

2 

where  N  is  the  buoyancy  frequency  defined  by 


(66) 


which  is  assumed  to  be  constant  throughout  the  stratosphere  and 
mesosphere.  The  term  RT'w/cpH  was  neglected  in  Lordi 1 s  formulation. 
Combining  (62)  and  (65)  as  a  single  expression,  we  can  write 


3T 1 

7 r 


=  r 


-  CONSTANT  x  /  D 


(67) 


i.e.,  temperature  can  be  expressed  as  a  group  of  terms  varying  in  time 
minus  terms,  constant  in  time,  times  an  integral  of  the  divergence. 

2.2.6  Continuity  Equation 

Expanding  the  horizontal  gradient  operator  of  equation  (47),  we 


can  express  the  continuity  equation  as  follows: 
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OP  /3t)  -  -  i  (V  •  VP  )do  -  J  P  Dda  -  *8PR 
*  0  *  0  *  H 


Multiplying  (68)  by  o/(oP  +p  )  yields 

*  R 


1  +  +  aP  1  w  p 

44  =  — -  /  (V  •  V  P  )do  -  *  [  /  Ddo  +  **  ]' 

Tt  oP  +p  o  *  aP  +p  q  nr 

*  R  *  R  * 


Using  (54),  (60),  and  (61),  we  can  isolate  divergence  and  mean  quanti¬ 
ties  to  obtain 

1  ♦  ,  w  p 

3<1  =  [  /  (V  •  VP  )do  +  P  6  R  R  ] 

Tt  aP  +p  o  *  *  M 

*  K 

1  _  w  p 

-  a0‘[  /  Dda  +  5  _R  R  ] 

0  H 

1  p  Z„  ZT 

-  o0  [  /  Dda  +  T  _R  exp(_R)  /  D  exp(:iZ)dz]  .  (70) 

o  H  H  z  H 

R 

The  continuity  equation  is  not  a  prognostic  equation  in  the  strato¬ 
sphere  and  is  used  primarily  to  compute  vertical  velocities  as  given 
by  (52).  Assuming  3q/3t  =  0  for  the  upper  part  of  the  model,  we  can 
still  write  a  combined  equation  such  as 

(3q/3t)  =  P  -  constant  x  /  D  (71) 

where  P  is  the  varying  part  and  described  by  the  first  two  terms  on 
the  right  of  (70)  for  the  troposphere  and  zero  for  the  stratosphere 
and  mesosphere.  Similarly,  the  constant  is  zero  for  the  upper  portion 
of  the  model . 


1 


2.2.7  Hydrostatic  Equation 

Rearranging  the  hydrostatic  equation  (8),  we  have 


3cP 


oP  +p 
*  R 


ii 

RT 
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(72) 


which  can  be  further  simplified  as  follows: 


3$ 

Uq  =  - 


RT  . 


(73) 


In  terms  of  layer  mean  quantities  and  deviations  from  layer  means,  we 
can  write 

3 *  =  -  RT  (74a) 

Iq 


3*1  =  -  RT'  •  (74b) 

W 

According  to  Kurihara  (1968)  the  form  of  (73)  yields  a  more  accurate 
estimation  of  geopotential  over  variable  topography  when  employing 
finite  difference  integrations. 

For  the  stratosphere,  we  have 

3*'  =  RT'  ,  (75) 

dF~  “FT 


with  a  similar  expression  for  layer  mean  quantities.  Equations  (74b) 
and  (75)  can  be  combined  into  a  single  expression  when  we  depict  the 
equations  in  the  discrete  state  later  in  this  chapter. 
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2 . 3  Vertical  Differencing 

In  this  model,  vertical  derivatives  and  integrals  are  computed 
via  finite  difference  methods.  Both  parts  of  the  model  employ  the  same 
differencing  techniques.  For  the  upper  portion  of  the  model  this  is  a 
slight  change  from  the  original  centered  space  differencing  technqiues 
as  used  by  Lordi  (1978). 

Since  vertical  velocity,  a  for  example,  is  formed  at  diagnostic 
levels  and  the  prognostic  variables  are  defined  at  alternate  levels, 
we  define  basic  differencing  schemes  involving  j  and  the  prognostic 
variables  denoted  by  X,  at  prognostic  level  r.  Similar  to  Corby  et 
al .  (1972),  these  relationships  are 

•  • 

(V  3<ii  =  X  r-b  ^  ;  (76) 

T7'r  r  - To - 


o  (X  -  X  )  +  a  (X  -X  ) 

.  3X  =  r+  b  r+1  r _ r-b  r  r-1  ; 

^ ’  To  ' r  2 Ao 


(77) 


o  (X  +X  )  -  o  (X  +-X  ) 

,3oX,  =  r+%  r+1  r  r- *s  r  r-1 

(ST"  V  - ZTo 


(78) 


Level  representation  and  indexing  are  denoted  in  Figure  1,  where 

Ao  =  .2  and  0  =  rAo.  When  r  M  we  use  the  following  relationship  at 

r 

the  boundary: 


[0  iU]D“  C*  d{ 

9o  K 


(79) 
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or  in  discrete  form 


,  4  z  [()-()] 

rw  1LJL-]  =  exp(  J<li}w  _  K  =  r=  \ 

LW  r\7  J  R  R  7  -  7 


(30) 


K  =  »s 


r=  <5 


Another  special  case  occurs  when  r  =  4, .  Here  o  =  0  and  (76  -  73) 

*  5 

can  be  simplified  accordingly. 

It  follows  from  the  specifications  of  (76  -  78),  that  the  finite 
difference  analog  maintains  the  following  derivative  relationship: 

3 ( ax )  =  *  3X  +  x  3^  •  (81 ) 

3  a  3  o  3  o 


It  can  also  easily  be  shown  that  a  similar  relationship  holds  for 


•  *2  *■ 

3(^V  )  =  y 

3o 


Op 

So 


*2 


y-  3_0 
dO 


(82) 


Vertical  integrals  are  treated  discretely  in  the  following 


manner: 


a  r 

(/  Xdo)  =  X  no  ,  (83) 

o  r+h  s=h  s 


where  r  represents  a  prognostic  level.  From  (83),  it  follows  that 
integral  quantities  are  formed  only  at  diagnostic  levels.  To  obtain 
values  at  the  prognostic  levels  centered  averaging  of  adjacent  diag¬ 
nostic  level  values  is  used. 

For  the  hydrostatic  equation  (74b)  for  the  troposphere,  we  employ 
the  following  discrete  scheme: 
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4>'  -  *'  =  -  Ra  T*  ,  (84) 

r+is  r-h  r  r 

where 

a  =  £n[(a  P  +p  )/  ( a  P  +p  )]  .  (85) 

r  r+*s  *  R  r-h  *  R 

Geopotentials  determined  by  (84)  correspond  to  diagnostic  levels  and 
centered  averaging  is  used  to  obtain  values  at  prognostic  levels. 

For  the  upper  portion  of  the  model,  (75)  is  expressed  discretely 
by  simple  centered  differencing,  i.e. 

RT 1 

-  4>  =  _ K  AZ  (86) 

K+  ij  K-  h  h 

and  averaging  is  used  to  obtain  prognostic  level  values. 

2.3.1  Mass  Conservation 

Lorenz  (1960)  showed  that  one  can  maintain  integral  constraints 
when  vertical  derivatives  are  replaced  by  finite  difference  analogs. 
This  was  the  reason  behind  defining  the  schemes  presented  in  the  pre¬ 
vious  section.  As  pointed  out  by  Arakawa  and  Lamb  (1977),  the  main¬ 
tenance  of  certain  integral  constraints  such  as  mass  and  total  energy 
conservation  by  finite  differences  may  not  be  critical  for  short  range 
numerical  weather  prediction  on  the  order  of  one  or  two  days.  However, 
long-term  integrations  beyond  a  few  days  require  conservation  of  these 
properties  to  prevent  non-linear  computational  instability,  false 
cascade  of  energy  into  smaller  scale  motions,  and  other  undesirable 
features  which  can  reduce  the  overall  forecast  accuracy  of  the  model. 
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In  this  section,  we  test  the  continuity  equation  for  the  main¬ 
tenance  of  mass  conservation.  Writing  (44)  in  discrete  form  as 
detailed  in  the  last  section,  we  have 


4%  >  +  w  p 

OP  /at)  =  -  H7  *  (p*v)ao)  -  JL1  .  (87) 

*  r=V  H 


If  we  integrate  (87)  over  the  entire  globe,  the  horizontal  gradient 
terms  vanish  and  mass  conservation  (aP^/at  =  0)  requires  that 

II  w  dS  =  0  ,  (88) 

S  R 


where  S  represents  the  global  surface.  Substitution  of  (54)  for  w 

R 

based  on  the  upper  model  formulation  clearly  indicates  that  (88) 
indeed  holds  and  shows  that  this  modeling  conserves  mass. 


2.3.2  Conservation  of  Total  Energy 

In  order  to  prove  total  energy  conservation,  we  need  to  derive 

the  governing  energy  equation  in  discrete  form.  To  derive  the  kinetic 

-► 

energy  equation,  we  multiply  the  momentum  equation  (5)  by  P  V  •  to 

*  r 

obtai n 


P  _ L  +  P  V  *V%v2  +  (P  V *o  iV)  :  -  P  v  -  P  RTV  *Vq  +P  V  *F.  (89) 

*  at  *  r  r  *  3ar  *rr*rr*rr 


-2 

Multiplying  the  continuity  equation  (7)  by  %V  yields 

r 


%v2(ap  /at)  +  %v2  v 


♦  %  3P  o 

(P  V  )  +  (%V2  *  )  =  0 

*  r  3a  r 


(90) 
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where  the  {  )  terms  in  (89)  and  (90)  involve  vertical  differencing 
r 

evaluated  at  prognostic  level  r.  However,  since  the  discrete  relation¬ 
ship  for  (82)  holds,  we  can  add  (89)  and  (90)  to  obtain  the  following 
kinetic  energy  equation: 


3P  *jV2  ^  +  +  3P  o%V2) 

*  r  +  v  •  (P  V  % V2 )  +  f  *  r  1 
3t  *  r  r  3o  r 


=  -  P  V  •  7*  -  P  RTV  •  Vq  +  P  V  •  F  .  (91) 

*r  r  *  r  r  *r  r 


In  a  similar  fashion,  we  can  multiply  the  thermodynamic  equation 

(6)  by  P  c  and  the  continuity  equation  (7)  by  c  T  and  add  the  result 
*  P  P  r 

ing  equation  using  the  discrete  property  of  (81)  to  yield 


3P  C  T  >  >  3£P  ac  T)  RT  P 

*  P  r  +  y  ■  (P  y  c  T  )  +  f  _  *.£J  =  _ Ll —  P  +  P  0  •  (92) 

at  *rpr  ao  r  o  P  +p  r  r 

r  *  R 


as  the  thermodynamic  energy  equation.  Adding  (91)  and  (92),  we  have  a 
total  energy  equation  of  the  following  form: 


3 [P  (%V2+C  T  )]  «.  +  +  3[P  o(%V2+c  T)] 

*  r  p  r  +  v  .  [p  V  (%V2+C  T  )]  +  {  _ L_  1 

at  *  r  r  p  r  9a  r 


(a) 


=  _  p  V  -V4>  - 

*  r  r 


(b) 


P  RT  V  -Vq  + 
*  r  r  r 


(c)  (d) 

RT  P  ♦  ♦ 

r  *  p  +  P  (V  *F  +Q  ) 
a  P  +p  r  *  r  r  r 
r  *  R 


(93) 


where  we  have  denoted  that  the  four  terms  on  the  right  of  (93)  by  (a) 
(b),(c),  and  (d)  for  future  reference. 
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Substituting  the  discrete  form  of  (50)  for  p  in  term  (c)  above 


allows  us  to  write 


(b)  +  (c)  =  -  P  RT  V  «Vq 
*  r  r  r 


p  py  ^  ^  w  p  ^  ^  r*  *  1  ^  ^ 

+  *  r  [a  V  *7P  -  R  R  -  i  (  l  7.(p  V  )Ao+  l  v*(p  v  )Ao)].  (94) 

a  P  +p  r  r  *  H  2  s=%  *  s  S-*  *  s 

r  *  R 


Since  the  first  term  on  the  right  of  (94)  cancels  with  the  first  part 
of  the  second  term  and  substituting  (43)  for  the  remaining  terms  in 
[  ]  yields 


P  RT  .  3P 

(b)  +  (c)  =  r  [P  a  +o  _J!J  (95) 

a  P  +p  *  r  r  at 
r  *  R 


Using  the  hydrostatic  equation  (8)  and  (95)  allows  us  to  write 


•  aP 

(b)  +  (c)  =  (-Poli)  +  (-  o  1*  __L)  .  (96) 

*  8or  So3t  r 


Substitution  of  (72)  into  the  second  term  on  the  right  of  (96)  gives 


us 


•  aP 

(b)  +  (c)  =  (-  P  a  i±)  +  (oRT  i!  __*)  .  (97) 

*  3o  r  3o  at  r 


The  horizontal  gradient  operator  of  term  (a)  can  be  expanded  as 
follows: 


(a)  =  -  v  •  (P  V  4>  )  +  $  7  •  ( P  V  )  . 

*  r  r  r  *  r 


(98) 


From  the  continuity  equation  (7),  the  above  term  can  be  written 
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-  ♦  3P 

(a)  =  -  7  •  (P  V  $  )  +  i>  (- _ _ 

*  r  r  r  3t 


P  i£) 

*  5o  r 


(99) 


or  with  further  manipulation  because  property  (81)  holds  for  the 
discrete  state,  (99)  becomes 


+  -  3P. 

(a)  =  -  v  •  (P  V  H  •  *  -JL 
*  r  r  r  at 


3P  oi> 

(  *  )  +  (p  a  ii)  . 


3o 


(100) 


Hence,  substitution  of  (97)  and  (100)  into  (93)  yields 


3 [P  (%V2  +c  T  )]  + 

*  r  p  r  +  v 

3t 


[P  V  (^V2  +c  T  )]  + 
*  r  r  p  r  r 


3[P  o(4V2+C  T+4)] 
{  _J1 _ P 


3  0 


I 

r 


3P 

=  -  (4>  -  oRT  i!)  — 1 
3a  r  3t 


+  P  (V  *F  +0  )  . 
*  r  r  r 


(101) 


If  we  integrate  the  first  term  on  the  right  of  (101)  from  the 
bottom  to  the  top  troposphere  boundary,  we  have 


4ir  3P 

Raj  3q_)  _Z  a o 
r=%  3o  r  3t 


3P 

-  _ *  Ao  I  ($  +$  +$  +$  +$  ] 

h  H  Zh  4% 


P_  (a  T  a  +o 
Aa  h  h  h 


T  a+aT  a+aT  a+oT  a 
1%  1%  1%  2%  2%  2%  3%  3*  3%  4%  4%  4% 


JJ  • 


(102) 


By  using  (84)  and  by  successive  substitution. 


we  can  write 


4  =  t  +  Ra  T  ; 

4  *  4%  4% 


* 

3 


4>  +  R(a  T  +a  T  ) 

*  4%  4%  3%  3% 


(103a) 


(103b) 
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4  =  4  +  R(a  T  +  a  T  +aT  ): 

2  *  4%  4%  3k  3%  2k  2k 


4  =  4  +  R{  a  T  +a  T  +a  T  +a  T  )  ; 
1  *  4%  4*  3%  3%  2%  2k  lk  lk 


♦  =  4  +  R(a  T  +a  T  +a  T  +a  T  +  a  T  )  ; 

B  *  4%  4%  3k  3k  2k  Ik  Ik  Ik  k  k 


(103c) 

(103d) 

( 103e) 


where  4^  =  4^_  the  surface  geopotential  and  4^  =  4  the  geopotential 
at  the  reference  pressure  boundary.  Using  (103a-e)  and  averaging 
diagnostic  level  values  to  obtain  prognostic  values,  we  have 

4  +4  +4  +4  +4 

k  Ik  Ik  3k  4  k 

=  54  +R(.5aT+1.5o  T  +2. 5a  T  +3. 5a  T  +4. 5a  T  )  .  (104) 

*  k  k  lk  lk  Ik  2k  3%  3k  4 k  4 k 

Substituting  (104)  into  (102)  and  using  the  actual  a  values,  (102) 
becomes 


4%  3P  3P 

l  -  (4  -  RoT  il)  _1  Ao  =  54  _JL  Ao  .  (105) 

r= k  3o  r  3t  *  3t 

Using  expression  (105),  we  can  now  write  the  entire  results  of  equation 
(101)  integrated  from  the  top  tropospheric  boundary  down  to  the  surface 
as 


,  4%  +  +  4  k  -*•„ 

°  [P  \  (%VZ+C  T  +4  )Ao  +  V  •  )  (^V^+c  T  +4  )Ao 

3t  *  r=k  r  p  r  *  r=%  r  p  r  r 


+ 


4 k  3[P  a(%V^+c  T+4) ] 

l  i  — _ E _ )  Ao 


r»% 


r 


4  k  * 

l  [P  (V  -F  +Q  ) ]Ao  .  (106) 

r=k  r  r  r 


da 


40 


Performing  the  discrete  differencing  as  specified  by  (78)  on  the 
vertical  flux  term  of  (106)  yields 


r  •  -2 

4%  3 [ P  o(%V  +C  T+$)3  . 

1  i _ *  . P _ j  Aa  «  -P  a  (%V^+c  T+$)n  (107) 

r=  %  3  a  r  *  R  p  R 


where  (  )  indicates  evaluation  at  the  reference  level.  If  we  substi- 
R 

tute  (107)  into  (106);  assume  adiabatic  and  frictionless  conditions; 
and  integrate  over  the  entire  global  surface,  we  obtain 


i_  //  [P  l  (*v2+c  T  +4  )Ao]dS  =//  p  0  (*v2+c  T+$)  dS  .  (108) 
3tS  r=%  rpr*  S*R  P  R 


A  similar  analysis  can  be  made  for  the  upper  portion  of  the  model 
to  yield 


P  25%  * 

—  If  [_§.  I  (%v2+c  T  +$  )exp(i£)Az]dS 
3t  S  H  K=%  K  p  K  R  H 


P  -*• 

=  //  J>  w*(%v2+c  T+4>)  dS  .  (109) 

S  H  R  P  R 


The  total  energy  for  the  entire  model  is  then  described  by  the  summa¬ 
tion  of  (108)  and  (109).  This  implies  that  total  energy  is  conserved 
only  if 

*  A  P  + 

//  P  0  (%V2+c  T+*)  dS  +  //  J>  w*(*V2+c  T+$)  dS  =  0  .  (110) 

s  *  R  p  R  s  h  R  p  R 


Since  it  follows  from  (4),  (42),  and  (53)  that  Pa  =  -  w*  s  , 

*  R  R  H 
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equation  (110)  does  indeed  hold  and  total  energy  is  conserved  through¬ 
out  the  entire  model. 


2.4  Semi -implicit  Time 
Differencing  Preliminaries 

Time  derivatives  are  predominantly  computed  via  semi-impl icit 
time  differencing.  This  method  slows  high  frequency  waves  that  would 
otherwise  limit  the  size  of  the  time  step  required  for  integration  of 
the  prognostic  equations.  Robert  (1969)  estimated  that  the  semi- 
impl  icit  time  method  can  use  a  time  step  6  times  that  of  a  fully 
explicit  time  differencing  scheme. 

Similar  to  Hoskins  and  Simmons  (1975),  we  can  write  the  finite 
difference  analogs  for  the  divergence  equation  (32),  the  thermodynamic 
equation  (67),  the  continuity  equation  (70)  and  the  hydrostatic  equa¬ 
tions  (74b)  and  (75)  in  matrix  form  as  follows: 


aCDD/at  =  [P]  -  v2([$']+R{TKq])  +  k  [v2d+2(LJ]  ; 

h  2 

a 


aCT' 3/at  =  [T]  -  t[d]  +  K  [v  (T'-T')]  ; 

h  e 


a C q 3 / a t  =  [P]  -  *[DJ  ; 


[*'-•']  =  Ry[T‘]  . 


(Ill) 


(112) 


(113) 


(114) 


where  [  ]  in  the  above  equations  Indicates  a  column  vector  of  length 
31  corresponding  to  the  prognostic  levels.  If  we  let  i  be  the  row 
index,  we  have  the  following  correspondence 


i 


i  =  1 


*  K  -  % 

■*  K  =  1% 

♦  K  =  2% 
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1  =  2 
1  =  3 


i  =  26  -  K  =  25% 

i  =  27  -  r  =  % 

1  =  28  *  r  =  H 


i  =  31  +  r  =  4%  . 

{T }  Is  a  31  x  31  matrix  composed  of  zeroes  everywhere  except  for 
values  on  the  lower  right  diagonal  positions,  x  is  a  31  x  31  matrix 
composed  of  terms,  constant  in  time,  which  multiply  divergence  to 
produce  temperature  tendency,  n  is  a  31  x  31  matrix  that  contains 
all  direct  effects  of  means  on  tropospheric  log-pressure  tendency. 

Y  is  a  31  x  31  matrix  composed  of  Az/H  and  a  terms  used  to  compute 
geopotential  hydrostatically.  Note  that  y  includes  the  averaging 
process  and  the  resulting  $'  values  of  (114)  are  prognostic  level 
values. 

If  we  let  j  be  a  summation  index  for  the  vector  terms  [P],  [T], 
and  [P],  which  corresponds  to  the  model  levels  as  does  i,  then  we 
can  define  these  more  complex  vectors  of  equations  (111  -  114)  using 
i  and  j  in  place  of  K  or  r  as  the  level  indicators.  For  example, 
i  =  28  really  means  r  =  H  and  so  on.  Hence,  we  have 
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Pi 


93  3A _  2 

1  [  i  -  cos  9  _ L  1  "  7 

T~  T5T  3$  1 


a  cos  9 


SB .  3A .  2 

1  „  [  1  -  cos  $  _1  1  -  v  E, 


"2  3X 
a  cos  9 


3  9 


(i  <  26) 


(115) 


(i  >26) 


T1  = 


-  1  r  3  ( UT  * ) i  +  cos  b  3  +  H, 

- 5“  L  5T~  3$  * 


(i  <  26) 


a  cos  9 


(116) 


-1  r  3(UT')i  +  cos  «  ]  +  H  .  M  >  26) 

2  3  X  39  r 

a  cos  9 


where 


RT 1  z . 

H  =  D  r  +  a(z  )(r-T')  -  — i-  !wu,+w*.jexp(r) 
K  i  i  i  e  i  2c  H  1+*  1  ’  H 

P 


expC^l)  ★  *  .  ,  . . 

*L[w.  ,(T:  ,-T! )  +  W.  (T  -T  )3  ’ 

2AZ  i+V  i  +  l  1  1  1-1 


ir .  y;  *  *  »V  jb  *  °i-*(VTi-in 


(VfVu  f  (v.-vp  )  -  i  (v  -?p,i*^_<o 

2p  1+*  j=27  0  *  j=27  J  Ao  1+*  H 


i  -*•  + 


P 
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(T.-T.  )  31  +  +  i-1  >  > 

— 2 — — -  ...  [  o  l  (V  ‘VP  )  -  l  (V  «7P  ) 
2P  *-*  j=27  J  *  j=27  J  * 


p*6' 

+  -  (a.  -1) 

A  a  1 


] 


RT!  1  >  +  i-1  +  +  +>wp 

- 1 _  [  *£  J,  (P  D  +V  *7P  }+^£  l  (P  D  +V  »7P  ) -o  V  *VP  +_Li] 

c  (o  P  +p  )  2  j  =27  *  J  J  *  2  j=27  *  J  J  *  1  1  *  H 

p  i  *  R 


RT  i  +  ^  i  - 1  wp 

-  -  1  [  —  }.  V  *7P  +  l  VVP  -  a  V  .VP  +  P  6'  R  R  ] 

c  (a  P  +p  )  2  j  =27  J  *  2  j=27  *  1  i  *  *  H 

p  i  *  R 


RT.  1  i-1  wp 

_ Le'  [^-1  0  +  l  D+T  R  R  ]  ; 

c  1  2  j =27  J  2  j=27  j  H 

P 


Next,  we  can  use  j  as  the  column  indicator  for  the  matrices  in  (111  - 
114),  corresponding  to  the  model  levels  similar  to  i,  to  define 


Mi  =  j  <  26) 


2 

t  =  (HN  Az/R)  + 

ij 

0  (elsewhere) 


expCiii-mz] 

H 

0  (elsewhere) 


+ 


is-expM  j-DAz] 
H 


(i=27 ,  j  <26  )  1 


0  (elsewhere) 


+ 


-T.) 

1 


2 


o  ( 31 > i >26 ,  j>26) 
i  +  % 

0  (elsewhere) 


1 (31>i>j>26) 

0  (elsewhere) 


6( o  -1 )  p 

0  (elsewhere) 


+ 


(T.-T.  1 

1  i-l 


a  ( i >27 ,  j>26) 
1-* 


0  (elsewhere) 


1  (i-l>j>26) 


0  (elsewhere) 


6  ( o  “1 )  p 

i-%  R  Az  expf-( j-l)Az-(Az/2)~| 


0  (elsewhere) 


(i  <  j  <  26) 


( 26<i <31 ,  j<*26) 


( i  >27 ,  j«;26) 

_i 
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RTi  - 
+  — L  e  ao 

c  1 

.5  ( 2 7 < 1 > j > 26 ) 

.5  ( i -1  >j  >26 ) 

.25  (i=j=27 ) 

+ 

i 

+ 

p 

0  (elsewhere) 

0  (elsewhere) 

0  (elsewhere) 

6  PR  Az  expf(-( j-UAz-Uz/2)n  (i>27,  j<26) 
TT~  Ao  H 


6  kR  Az  exnfM  j-l)Az-(Az/2)i  (i=27,  j<26) 
T  H  Ao  T1 

0  (elsewhere) 


(118) 


77 

ij 


0  (i<26,  j)  + 


o  3  Ao  ( i >26 ,  j>26) 
i  i 


0  (elsewhere) 


_ p 

a  3  6  R  AZ  exnr-(j-l)A2-(Az/2)l  (i>26,  j<26) 
i  i  IT-  FI 

0  (elsewhere) 


(119) 


"  .5  ( i=j<26) 

x 

1  ( i > j  <26 ) 

T 

_0  (elsewhere) 

0  (elsewhere) 

Y  .  =  “ 
10  H 


a  ( i  «s26 ,  j >26) 

.5a  (i=j>26) 

a  (i<j>26) 

j 

j 

j 

+ 

+ 

+ 

0  (elsewhere) 

0  (elsewhere) 

0  (elsewhere) 

(120) 


.J 
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Using  the  semi -imp! icit  scheme  of  Robert  et  al.,  (1972),  we 
define  the  following  notation: 


t+At  t-At 

<5  X  i  x  'x  ,  (121) 

t  2  At 


t  t+At  t-At 

X  H  x  ~x _  ,  (122) 

2 


where  the  superscript  represents  the  time  step.  Using  (121)  and  (122), 
we  can  express  prognostic  expression  (111  -  113)  in  semi-implicit  form 
as  follows: 


6  [D]  =  [p]  -  +  R  { T  jCq]1)  +  K  [V2D+2(D/a2  )]t_At  ;  (123) 


t  2  t-At 
<5  [T1]  =  [r]  -  x[D]  +  K  [V  (T'-T1)]  ; 
t  he 


(124) 


$  [q]  =  [P]  -  xi[D]  . 

Since  we  can  rearrange  (121 )  and  (122)  to  obtain 

t  t-At 


X 


6 


X 

t 


=  X  +  At  <5  X 

t 


t  t-At 
=  X  -X 

At 


(125) 


(126) 


(127) 


we  can  substitute  using  (126)  and  (127)  to  express  (123)  as  follows: 
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.t  ,  t-At 


[oi'-cor  .  m-  vz([*']t'‘'t+  out[»,j*RiT|[<i]t"“*itl,iTi5t[<,]) 

At 


2  2  t-At 

+  K  [V  D  +  2 (0/a  H 
h 


(128) 


Before  we  can  proceed  further,  we  need  to  develop  an  expression 

the  continuous  time  derivative  of  (114), 


for  6  [*'  ].  First,  we  take 
t 

which  yields 


(129) 


30']  =  at*;]  +  Ry  3tT'  3  +  R  IT’]  • 
at  at 

I„  order  to  compute  »,/«,  we  need  to  be  .bit  to  compute  t.  From 

the  definition  of  a  given  by  (85), we  can  write 


8a 


I  „  L_  [into  P  +P.)  -  »nlo  P  +PR^  ‘ 


(130) 


at  at 


r+h  *  R 


r-% 


Differentiating  the  right  side  of  (130)  yields 


3o 


at 


r+% 


r-% 


ap , 


]  TT-  * 


(131) 


o  p  +p  'a  P  +P  3t 
r+k  *  R  r-%  *  R 


Substitutions  of  the  discrete  form  of  («>  for  *V«  «»  «31>-  •”«» 
us  to  write 

4% 

o  J  ,  -*rPr 

(  r»»>  _Iil_  ^  -JLS-  _  2. 


8a 


at 


V-V  ^ 


(P  D  +V  *VP  ) Ao  ]  .  (132) 

*  s  S  * 


A 


Since  all  the  other  terms  of  y  are  either  zero  or  constant  in  time,  the 

only  remaining  terms  of  3y/3t  involve  9a  /9t  which  can  be  expressed 

r 

diagnostically  by  (132). 

If  we  define 


r  ;  3y 

at 


(133) 


Y  =  Y(o)  +  y'  , 


(134) 


we  can  express  (129)  as  follows: 


ai‘i|1]  +  Ry'  3[t‘3  +  Ry  3^t> ]  +  Rr  [t1]  • 

3t  at  at  at 


(135) 


If  we  compute  3[T']/at  for  the  second  term  on  the  right  of  (135)  by 
computing  the  right  side  of  (112)  directly  and  call  this  value  [e], 
then  we  can  express  (135)  in  discrete  form  as  follows: 


6  [*']  =  6  O']  +  Ry'[8]  +  Ry  5  [T 1  ]  +  RF[T '  ]  .  (136) 

t  t  *  t 


Substituting  (124)  for  6  [T1]  in  (136)  and  substituting  (136)  for 
6  O']  in  (128)  and  substituting  (125)  in  (128)  yields 


{  I  -  RAt2(7x  +  (T}ir)V2}[D]t 


[0]t_At  +  AttfP]  -  V2([$']t'At  +  R{T }[q]t_At )  } 


-  At2V2{6^[<t;]+R(Y,[0]+r[T,]+7lT]+{Tj[F]+YKh[V2(T,-T^)]t"At  j 
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2  2  t-At 

+  At  K  [V  D  +  2(D/a  )] 
h 


(137) 


where  I  is  the  identity  matrix.  With  the  spectral  representation  to  be 

introduced  in  the  next  section,  we  will  be  able  to  solve  (137)  for  [ D ]  t 

t+At  t+At  t'+At 

and  subsequently  solve  for  [D]  ,  [T * ]  ,  and  [q] 


2.5  Spherical  Harmonic  Representation 
Except  for  vertical  derivatives  and  the  computation  of  non-linear 
terms,  which  are  formed  at  grid  points,  other  model  computations  are 
handled  spectrally.  The  variables  x>  *'»  U,  V,  T1  and  q  can  be 
expanded  in  terms  of  spherical  harmonics  as  follows: 


2  +J  lrol+L  m  m  m,  m 

l*»  x,  *  }  -  a  1  l  {*  ,  x  .  *  }Y  *. 

m=-J  £=  1ml  *■*■*■*■ 


+J  lml+L+1 


{U.  V}  =  a  X  l  {U1;,  V"}y"  ; 

m=-J  £=|m|  4  1  l 


+J  Iml+L 


m,  m 


(PJ  =a  l  l  {P  }Y  ; 
m=-J  A= 1 m I  *  * 


+J  lml+t  m  m.  m 

iT'.qj  =  l  l  {T  .  q  iY  , 
m=-J  t=|ml  *  *  * 


where 


m  m  imx 
Y  =  P  (u)e 

l  i 


(138) 


(139) 


(140) 


(141) 


(142) 


v  -  sin  <j>  . 


(143) 
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m 

P  (p )  is  an  associated  Legendre  polynomial  normalized  to  unity  with  m, 
l 

the  longitudinal  wavenumber,  and  i,  the  latitudinal  index.  L  deter¬ 
mines  the  limit  of  the  parallelogramic  truncation.  The  {  j  terms  on 
the  right  of  (138  -  141)  represent  the  respective  spherical  harmonic 
coefficients.  J  is  the  longitudinal  wavenumber  truncation. 

After  the  vertical  derivatives  and  non-linear  products  have  been 
formed  at  grid  points,  the  resulting  terms  can  be  transformed  in  terms 
of  Fourier  series  as  follows: 

,  ,  imx 

(Z,  V,  T,  P,  Y ' C 8] ,  r[T' ]}  =  l  {Z  ,  V  ,  T  ,  P  ,  6  ,  K  }e  (144) 

m=_j  m  m  m  m  m  m 


where  the  m  subscript  represents  the  respective  Fourier  coefficient. 
The  Legendre  transform  defined  by 

m  7r/2  m 

(  )  =  /  (  )  P  (p)cos*  d*  (145) 

l  _t,/2  m  i 


can  be  applied  to  the  Fourier  coefficients  in  (144)  to  obtain  the 
appropriate  spherical  harmonic  coefficients.  Since  in  the  model,  the 
Fourier  series  defined  by  (144)  are  formed  at  Gaussian  latitudes, 

(145)  can  be  computed  exactly  for  each  term  up  to  the  point  of 
truncation. 

One  of  the  most  important  advantages  in  using  spherical  harmonics 
is  that  horizontal  derivatives  can  be  computed  accurately  and  easily. 
The  following  derivative  properties  result  from  using  spherical 
harmonics: 


(146) 
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truncation,  a  significant  advantage  over  grid  point  methods. 

Using  the  spectral  relationships  outlined  so  far  in  this  section, 
we  can  write  (27),  (111  -  113)  as  the  following  spectral  equations: 
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(157) 


while  the  [  ]  represents  a  column  vector  as  defined  in  the  previous 
section.  Note  that  the  hydrostatic  equation  is  treated  discretely  and 
geopotential  values  computed  from  (114)  are  then  transformed  for  input 
to  the  spectral  divergence  equation. 

2.6  Conclusion  of  Semi-implicit  Differencing 

In  spectral  form  using  relationships  developed  in  the  previous 
section,  we  can  express  (137)  as  follows: 
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According  to  Hoskins  and  Simmons  (1975),  an  equation  of  the  form  of 
(159)  is  essentially  a  gravity  wave  equation.  We  may  invert  the 
matri x 


2 

M  i  {I  +  RAt  aU+1)(7t+{T}h)}  ,  (160) 

2 

a 

since  (7t+{T)")  has  positive  eigenvalues.  Since  M  is  constant  in  time, 
-1 

M  needs  only  to  be  calculated  once  before  the  start  of  prognostic 

_  mt 

integrations.  Hence,  we  can  solve  for  [D£ J  on  the  left  side  of 
(159). 

m  t  m  t-At 

Using  this  value  of  [D^]  ,  we  can  compute  [D^J  from  (122), 

[Tm]t+At  and  [qm]t+At  from  the  spectral  equivalents  of  (124)  and  (125), 
l  i 

i .e. ,  we  have 
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Since  we  cannot  include  the  vorticity  equation  in  the  semi- 
implicit  scheme  since  it  is  not  related  to  divergence,  we  use 
centered  time  differencing  of  (154)  with  direct  solution  as  follows: 
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Equations  (161  -  164)  yield  prognostic  values  for  D  ,  T  ,  q  ,  and  S 

X  X  X  X 

from  which  we  can  compute  new  values  of  ^  and  X™  from  (152  -  153)  and 

In  turn  compute  Um  and  Vm  from  expressions  (150  -  151).  We,  then,  can 
i  i 

construct  real  values  for  0,  c,  T‘,  U,  V,  and  q  from  the  spectral 
coefficients  and  subsequently  compute  p+  and  hydrostatically  compute 
*. 


56 


2.7  Modified  Euler  Backward  Scheme 
In  order  to  start  model  integrations  from  an  initial  state,  we  use 
the  modified  Euler  backward  time  differencing  scheme  after  Kurihara 
(1965).  If  we  let  X  represent  any  prognostic  variable  and  Y  represent 
the  right  side  of  the  prognostic  equation,  i.e., 

3X/3t  «  Y  ,  (165) 

we  can  represent  the  three  steps  of  this  iterative  process  as  follows; 

At. 


( t+  7")  t  At  t 

(166) 

Step  1  : 

X  =  X  +  7-  Y  , 

Step  2: 

At 

(t+At)*  t  (t+£“) 

X  =  X  +  At  Y 

* 

(167) 

Step  3: 

(t+At)  t  (t+At)* 

X  =  X  +  At  Y 

(168) 

This  method  of  time  differencing  is  applied  directly  to  the  spectral 
equations  (154  -  158).  Since  the  scheme  is  a  damping  one  horizontal 
diffusion  terms  are  neglected  when  it  is  used.  The  modified  Euler 
backward  scheme  also  removes  any  computational  modes  generated  from 
centered  time  differencing  (Haltiner,  1971)  and  consequently  is 
interspersed  every  37  time  steps  in  place  of  the  semi-implicit  scheme. 


CHAPTER  3 


SPECIFICATION 

3.1  Model  Initialization 

Initial  fields  of  geopotential,  temperature,  stream  function  and 
log-pressure  must  be  specified  before  the  start  of  model  integrations. 
For  these  initial  conditions,  we  specify  zonal  profiles  that  are  non- 
divergent  and  essentially  have  a  non-di vergent  tendency.  These  latter 
conditions  are  essential  so  that  high  frequency  oscillations,  which 
could  adversely  affect  the  model's  behavior,  are  basically  eliminated. 
In  this  section,  we  will  describe  the  initialization  process  used  to 
produce  the  initial  profiles  used  in  this  model. 

Our  goal  initially  was  to  use  the  initial  conditions  of  Lordi 
(1978)  for  the  upper  portion  of  the  model,  since  they  met  the  criteria 
specified  above.  However,  any  realistic  tropospheric  fields,  built  to 
match  these  conditions,  resulted  in  convectively  unstable  conditions  at 
the  interface  between  the  models  near  the  equator  and  some  alterations 
had  to  be  made.  By  subsequently  adjusting  the  initial  temperature 
deviations  in  the  lowest  three  levels  of  the  stratosphere  south  of 
45N,  we  were  able  to  obtain  a  temperature  field  that  was  convectively 
stable  and  also  cl imatologically  more  agreeable. 

To  arrive  at  a  balanced  state  in  the  troposphere,  we  initially 


started  with  constant  pressure  level  zonal  temperature  deviations 
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from  level  means,  based  on  January  climatic  tables  from  Oort  and 
Rasmusson  (1972).  Data  for  polar  latitudes  were  extrapolated  from 
January  1974  and  1976  data  from  National  Meteorological  Center  (NMC) 
observational  grids.  The  above  temperature  fields,  which  were 
specified  at  5°  latitude  increments,  were  first  linearly  interpolated 
to  Gaussian  latitudes  used  in  the  model.  Using  Lordi's  original  lower 
boundary  geopotential  deviation  values  as  an  upper  boundary  condition, 
we  integrated  hydrostatically  downward  to  obtain  geopotential  at 
intermediate  pressure  levels  to  a  pre-selected  constant  pressure 
surface,  denoted  ]j  ,  which  corresponds  to  the  mean  surface  pressure  at 
our  o  =  1  level.  By  assuming  T  to  be  constant  (  =  T  )  with  height  near 

this  surface,  we  can  integrate  the  hydrostatic  equation  from  p^,  the 

actual  surface  pressure,  to  and  obtain 

*n(p  /p)  1  (*-*-)  .  (169) 

*  *  RT 

* 

Now  let  tn  =  *1)  +  4>'d  and  *0  =  Tp  +  <i>'7T  ,  where  the  bar  denotes  a 
layer,  mean,  and  assuming  that  ,  and  $*D  =  0,  corresponding 

to  a  flat  surface,  we  obtain 

p  =  p  exp($'—  /RT*)  ,  (170) 

*  *  M* 

where  temperature  values  were  determined  by  interpolation. 

Using  the  resulting  p+,  we  can  compute  the  pressure  corresponding 
to  the  five  tropospheric  prognostic  levels.  Next,  we  interpolate  the 
constant  pressure  level  temperatures  to  these  sigma  levels  and  inte¬ 
grate  hydrostatically  downward  in  the  sigma  system  from  the  interface 
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boundary  geopotential  used  previously  to  obtain  geopotentials  at  diag¬ 
nostic  levels  of  the  troposphere.  However,  because  of  interpolation, 
the  surface  geopotential  that  resulted  from  this  integration  was  not 
identically  zero  as  desired.  We  force  this  solution  by  redistributing 
the  differential  at  the  surface  logarithmically  to  higher  levels  and 
adjust  the  temperatures  at  prognostic  levels  to  maintain  hydrostatic 
balance.  These  temperatures  now  allow  us  to  integrate  the  hydrostatic 
equation  upwards  as  is  done  in  the  model  and  obtain  boundary 
geopotentials  identical  to  Lordi's. 

Averaging  the  diagnostic  level  geopotentials  to  obtain  prognostic 
level  values  as  is  done  in  the  model  and  using  the  corresponding 
temperature  and  pressure  fields,  we  can  then  quadratical ly  solve  the 
gradient  wind  balance  equation  for  our  o-coordinate  system 

uff  +  u  tan  .  1  3*'  .  0R(T  +  T 1 )  1  3P*  .  (171) 

a  a  3$“  <jP^  +  pr  a  3$ 

for  the  zonal  velocity  field  u.  The  mean  temperature  profile  T  is 
assumed  to  be  that  of  the  standard  atmosphere  at  45°N. 

Since  reflective  symmetry  across  the  equator  is  assumed,  the 
resulting  zonal  wind  field  may  be  expanded  at  prognostic  levels  in 
terms  of  the  following  series: 

u($,oJ  =  l  Ar(o)cos  rp  (172) 

r=l,R,2 


where  R  =  33  and 


Ar(o) 


4  "/2 

Z  J  uU,o)cos  r$  d<4>  , 

IT  o 


(173) 
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satisfying  the  boundary  conditions  that 

u(t  n/2,  o)  -  0  .  074) 

Integration  of  (173)  was  performed  numerically  using  a  Simpson's  rule 
adaption  for  unequally  spaced  arguments. 

Since  horizontal  divergence  is  assumed  to  be  zero,  it  follows 
from  the  continuity  equation  that 

(3y/3y)  =  -  a  l  Ar(o)cos  ry  ,  (175) 

r=l,R,2 

which  can  be  integrated  assuming  a  boundary  condition  that  y ( 0 , o )  =  0 
to  yield 

*U.°)  »  -  a  l  Ar(o)(l1n...ri.)  (176) 

r=l ,R,2  r 


We  can  then  transform  the  resulting  y  from  (176)  to  obtain  initial 
m 

y  fields  for  the  troposphere. 
i 

Following  the  model  procedures,  we  can  determine  U  and  ;  from 

m  m  m  m 

y  .  We  also  compute  the  spectral  coefficients  *  ,  T  ,  and  q  for 
l  l  i  i 

geopotential,  temperature,  and  log-pressure,  respectively.  Assuming 

m 

D  is  zero  and  neglecting  diffusion  and  friction,  the  divergence 
i 

tendency  equation  (156)  can  be  solved  to  yield  spectral  divergence 

-10  -15 


tendencies  on  the  order  of  10 


to  10 


In  order  to  reduce  this 


m  ,  4  m 

tendency  even  further,  we  set  3D  /3t  to  zero  and  solve  (156)  for  y  . 

i  l 

We  then  compute  *  from  the  coefficients  and  compute  T'  hydrostatically 

after  reversing  the  averaging  process  and  obtaining  diagnostic  level 

m 

geopotential  fields.  Transforming  T'  to  and  back  to  T',  we  then 


J 
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recompute  £  at  diagnostic  levels,  average  £  to  obtain  prognostic  level 

m 

values  and  then  compute  £^.  We  then  solve  (166)  for  divergence 

tendency  and  again  compute  a  £m  to  make  this  tendency  zero.  We  then 

i 

continued  with  the  remainder  of  the  process  outlined  above  and  itera¬ 
tively  continued  this  process  six  times  until  apparent  convergence  was 

-14  -18 

reached,  yielding  divergence  tendencies  on  the  order  of  10  to  10 
This  was  considered  sufficiently  small  so  that  any  initial  imbalance 
would  not  unduly  influence  model  integrations. 

As  mentioned  earlier,  this  initialization  procedure  resulted  in 
convectively  unstable  conditions  at  the  interface.  The  boundary  level 
geopotential  also  no  longer  exactly  matched  that  of  Lordi's  because  of 
the  above  balancing  process.  Both  conditions  mandated  changes  from 
Lordi's  initial  conditions.  To  initialize  the  upper  model  levels,  we 
started  with  Lordi's  initial  temperature  distibution  and  restructured 
the  temperatures  at  the  lowest  three  levels  at  lower  latitudes  in  keep¬ 
ing  with  observational  data  and  this  resulted  in  a  convectively  stable 
situation. 

Using  the  new  temperature  profile,  we  then  integrated  hydrostati¬ 
cally  as  in  the  model  using  the  boundary  geopotential  that  results 
from  the  tropospheric  formulation  to  obtain  diagnostic  level  geopoten¬ 
tials  which  were  averaged  to  prognostic  levels.  Using  the  gradient 
balance  equation  for  the  upper  portion  of  the  model,  we  solve 

u(f  +  u  tan  =  _  1  3£' 

a  a  If  (177) 


quadratically  for  u  and  can  proceed  as  with  the  troposphere  in 
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obtaining  an  initial  stream  function  field.  We  then  can  solve  (156) 

for  the  stratosphere  and  obtain  divergence  tendencies  on  the  order 

-11  -14  m 

10  to  10  .  Computing  a  residual  4>^  as  with  the  troposphere,  we 

can  find  a  corresponding  T™  field  from  the  hydrostatic  relationship 

-25 

yielding  tendencies  on  the  order  of  10  or  nearly  zero.  This  bal¬ 
ance  is  more  effective  in  the  stratosphere,  since  (177)  involves  only 
a  balance  between  u  and  <p,  whereas  (171)  also  involves  P*,  T  and  T'. 

The  above  analysis  yielded  initial  fields  for  4  ,  T  ,  0/  and  q 

l  l  l  i 

that  are  essentially  non-di vergent  and  possess  negligible  divergence 
tendency.  The  zonal  wind  and  temperature  profiles  for  this  case, 
which  we  will  denote  as  Case  1  are  depicted  in  Figures  2a  and  3a, 
respectively.  The  upper  part  of  these  profiles  matches  Lordi ' s  to  a 
high  degree. 

We  also  wanted  to  examine  the  case,  denoted,  Case  2,  which  had  a 
much  stronger  polar  vortex  in  the  stratosphere  and  mesosphere.  Using 
the  same  troposphere  as  in  Case  1,  we  adjusted  the  upper  model  tem¬ 
peratures  by  the  following  factor: 


T'  =  T'  [1 

CASE  2  CASE  1 


+  (1 


exp{-(J^]/i)  j]  , 


(178) 


where  K  represents  the  stratospheric  level  index.  This  formulation 
essentially  doubles  temperature  deviations  at  levels,  six  and  higher, 
above  the  interface  and  increases  more  gradually  the  intermediate 
levels  just  above  the  interface.  Using  these  temperatures,  we  can  then 
initialize  for  Case  2  using  exactly  the  same  upper  model  initializa¬ 
tion  as  Case  1.  Initial  zonal  wind  and  temperature  fields  for  Case  2 


L. , 


63 


are  depicted  in  Figures  2b  and  3b. 

Mean  layer  values  for  geopotential,  "fc,  although  not  required  in 
stratospheric  computations,  are  required  in  the  troposphere  to  compute 
height  fields  used  in  vertical  differencing  as  in  (80);  to  compute 
temperature  lapse  rates  for  use  in  convective  adjustment  procedures 
to  be  discussed  later  in  this  chapter;  and  to  interpolate  o-level 
fields  to  constant  height  levels  for  display  purposes.  These  values 
were  computed  by  using  Tp  as  a  bottom  boundary  and  standard  atmosphere 
temperatures  for  45N,  integrating  similar  to  (84)  and  center  averaging 
the  result  to  the  prognostic  levels. 

3.2  Orographic  Forcing 

Surface  geopotential  must  be  specified  as  a  lower  boundary  condi¬ 
tion  in  the  model.  It  is  required  by  the  hydrostatic  equation  (114) 
and  the  semi -implicit  divergence  equation  (159).  As  mentioned  in  the 
last  section,  initial  conditions  were  based  on  the  assumption  of  a 
flat  surface  at  the  lower  boundary.  We  tested  the  model  and  the 
balanced  state  of  these  conditions  by  keeping  the  surface  unchanged; 
neglecting  friction,  diffusion,  and  heating/cooling;  and  integrating 
the  model  for  10  days.  During  this  period,  divergence  increased  only 
slightly  and  the  initial  zonal  fields  remained  virtually  unchanged, 
indicating  the  equilibrium  state  of  the  initial  fields  and  the 
capability  of  the  model  to  maintain  them  in  the  absence  of  external 
forcing.  However,  our  interests  are  on  the  effects  of  orographic 
forcing  on  the  troposphere  and  stratosphere.  In  order  to  study  these 
effects,  we  slowly  introduce  terrain  features  at  the  lower  boundary 


which  are  reflected  in  geopotential  deviations  at  the  surface  that 
force  the  model  levels  above. 
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We  assumed  a  very  idealistic  and  simplistic  topographic  pattern 
to  simulate  the  continental  land  mass  distibution  of  the  Northern 
Hemisphere.  Due  to  spectral  truncation  used  in  the  model,  only  the 
zonal  ultra-long  waves  1  and  2  were  included.  The  longitudinal 
distribution  of  the  surface  height,  z+,  at  45N  can  be  expressed 
empirically  as  follows: 


z  (45N)  =  A  sin(m  X)  +  A  sin(m  X  -  7I)  +  z  ,  (179) 

*  1  1  2  2  7  * 


where  A  and  A  are  the  amplitudes  of  zonal  wave  numbers  m  =  1  and 
1  2  1 

m  =  2,  respectively;  7  is  the  initial  flat  surface  height  and  was 
2  * 

used  to  determine  the  initial  mean  surface  pressure  p*  based  on 
standard  atmosphere  profiles.  The  latitudinal  variation  can  be 
expressed  by 


z 

★ 


(*)  = 


z*(45N)sin 


(2$) 


(180) 


For  the  model  runs  to  be  described  in  this  thesis,  we  assumed 

7  =  A  =  240  m  and  A  =  150  m.  Although  there  is  no  mathematical 
★  2  1 

requirement  for  z^>  0,  the  above  choice  of  7^  insures  this  condition. 

The  45N  longitudinal  representation  of  the  z^  field  is  depicted  in 

Figure  4.  The  maximum  z#  for  the  simplistic  European/Asian  land  mass 

is  630  m  (A  +A  +7  )  and  for  the  North  American  continent,  it  is 
1  2  * 

slightly  more  than  half  that  at  330  m  (A  -A  +7  ).  Although  these 

12* 

magnitudes  are  small,  when  compared  to  those  used  by  Bourke  et  al . 


65 


(1977)  in  their  spectral  model,  it  should  be  remembered  that  their 
model  included  much  higher  harmonics  and  their  spectral  representations 
for  topography  were  able  to  be  more  confined  than  in  our  model. 

However,  the  relative  magnitude  differential  between  the  peaks  for  the 
continents  is  comparable.  Model  runs  indicate  that  the  forcing  used 
is  sufficient  to  generate  geopotential  forcing  at  the  tropopause  on 
the  order  of  300  gpm  by  combining  the  wave  1  and  2  components.  This 
easily  exceeds  the  wave  amplitude  of  200  m  at  300  mb,  which  Schoeberl 
and  Strobel  (1980a)  say  are  needed  for  stratospheric  warmings  to  occur. 

In  order  not  to  shock  the  model  by  the  sudden  introduction  of 
topography,  it  was  slowly  turned  on  according  to  the  following 
al gori thm 


=  g(z^  -  zj[1  -  expf-t/t^)]  ,  (181) 

5 

for  time  t  and  t  =  2.5  x  10  s.  With  this  formulation,  4  attains 
o  * 

about  90%  of  its  maximum  value  after  7  days.  Since  i>^  is  therefore  a 
function  of  time,  we  can  obtain  the  time  dependent  values  required  by 
(159)  and  the  current  time  values  required  by  (114). 

3.3  Convective  Adjustment 

Motions  generated  by  models  become  gravitationally  unstable  if 
the  temperature  lapse  rate  exceeds  that  of  the  dry  adiabat  given  by 

y  =  9/c  .  (182) 

d  p 


To  prevent  this  instability,  we  adjust  forecast  temperatures  when  they 


exceed  this  critical  value.  This  process  is  known  as  convective 
adjustment.  8ecause  of  the  inherent  approximation  in  spectral  models 
of  going  from  real  space  to  transform  modes,  some  accuracy  is  lost  to 
truncation.  For  this  reason,  we  define  a  slightly  weaker  lapse  rate 
as  our  adjustment  criteria,  namely, 

y  =  -  ( g-.5)/c  .  (183) 

s  P 


This  allows  for  roundoff  errors  of  the  transform  process. 

In  the  model,  the  lapse  rate  is  first  checked  at  the  interface. 

If  (183)  is  exceeded,  an  adjustment  is  performed.  Levels  in  the  upper 
atmosphere  are  then  checked  and,  if  necessary,  adjusted  sequentially 
upward.  The  remaining  tropospheric  layers  are  then  processed  downward. 
This  adjustment  process  consists  of  the  following  steps: 

(1)  Compute  lapse  rate  between  adjacent  prognostic  levels,  i.e. 


(T  +  T')  -  (T  +  T 1  ) . 

AT  =  i+l  1  . 

AZ  Z  -  Z 

i+l  i 


(184) 


As  elsewhere  in  the  model,  T  values,  based  on  the  standard  atmosphere 
at  45N,  are  used  for  layer  means  based  on  the  mean  surface  pressure 
of  the  layer  and  are  initially  specified  before  the  start  of  model 
integrations.  For  the  stratosphere,  we  use  the  z  values  given  by  (1). 
For  the  troposphere,  we  compute  height  by  the  following: 


z  =(*+*' )/g  ,  (185) 

r  r 


for  any  prognostic  level  r. 
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(2)  If  y  is  not  exceeded,  we  check  the  remaining  levels.  If  y 
s  s 

Is  exceeded,  we  must  insure  that  any  adjustment  made  on  temperature 
preserves  the  internal  energy  between  the  levels  being  adjusted.  This 
means  that 


P 

s  exp{-  £)Tdz  =  constant  , 


(186) 


/  c  P  Ido  =  constant  , 
P  * 


(187) 


for  the  upper  and  lower  portions  of  the  model,  respectively.  The 
constraints  denoted  by  (186)  and  (187)  imply  the  following: 


p  -z  P  -z 

.  s  exp(  *~^)T*  az+P  T*  Ao  =  s  exp(  K=J*)T  az+P  T  Ao  , 
T r  H  K=%  *  r»%  FT  ^  K=%  *  r=% 


(188) 


K+1.)t* 

K  K+l 


"  eXp(T)TK  " 


-z 

+  exp(  k)T 


(189) 


T*  +  T*  =  T  +  T  ,  (190) 

r+1  r  r+1  r 

for  the  interface  levels,  upper  model  levels,  and  tropospheric  levels, 

respectively,  where  T*  represents  the  adjusted  temperature  and  T  the 

original  temperature  at  the  level  noted  by  the  subscript. 

(3)  If  adjustment  is  required,  we  specify  that  the  new  lapse  rate 

equals  y  .  For  example, 
s 

T*  -  T*  «  y  (z  -z  )  (191) 

r+1  r  s  r+1  r 


i 
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for  the  troposphere.  Combining  (190)  and  (191),  we  can  determine  T* 

r+1 

and  T*.  Similar  determinations  can  be  made  at  the  interface  and  for 
r 

the  upper  model  levels. 

(4)  Finally,  after  all  levels  are  checked,  if  an  adjustment  were 

necessary,  geopotential  is  recomputed  using  (114)  and  spectral  coeffi- 
m  m 

dents  for  T  and  4  are  computed  based  on  the  new  values.  The  above 

i  i 

process  insures  that  convective  stability  is  maintained  and  that 
thermal  energy  is  conserved. 

3.4  Spectral  Truncation  and 
Horizontal  Diffusion 

Following  Lordi  (1978),  parallelogramic  truncation  is  used  as 
indicated  In  (138  -  141)  with  J  =  4  and  L  =  24.  We  also  assume  sym¬ 
metry  across  the  equator,  which  means  that  4,  C,  Z,  and  V  are  expanded 
in  asymmetric  harmonics  and  the  other  spectral  variables  need  only  be 
expanded  in  symmetric  spherical  harmonics.  This  cuts  in  half 
computer  core  storage  that  would  otherwise  be  needed  to  store  these 
spectral  coefficients. 

The  above  truncation,  along  with  the  fact  that  in  real  space 
there  are  N  =  18  points  around  each  latitude  circle  and  I  =  20 
Gaussian  latitudes,  insures  both  alias-free  calculations  of  the 
Fourier  coefficients  by  (144),  which  makes  the  transform  procedure 
for  evaluating  non-linear  products  equivalent  to  direct  spectral 
multiplication,  and  also  integrals  with  respect  to  latitude  may  be 
computed  exactly  by  means  of  Gaussian  quadrature.  For  details,  see 
Lordi  (1978). 
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As  mentioned  In  Chapter  2,  horizontal  diffusion  was  introduced  to 
inhibit  the  spectral  blocking  that  would  exist  because  of  the  extensive 
horizontal  truncation  used  in  this  model.  We  specified  a  horizontal 
diffusion  coefficient  as  follows: 

6  2  -1 

K  =1.5x10ms  ,  £<17 

h 

6  2-1 

K  =  6.4  x  10  m  s  ,  £  >  17  . 

h 

However,  diffusion  was  also  only  selectively  applied  to  the  higher 
harmonics  meeting  the  following  criteria: 

m  =  0  ,  £  >  17 

m  =  1  ,  £  >  15 

m  =  2  ,  £  >  13 

m  =  3  ,  £  >  11 

m  =  4  ,  £  >  9  . 

For  the  levels  at  and  above  43.5  km,  a  fourth  order  filter  developed 
by  Shapiro  (1971)  is  also  applied  to  eliminate  aliasing  that  can  exist 
in  the  vertical  due  to  finite  differencing.  Because  the  small  vertical 
modes  propagate  rapidly  upward  and  amplify  in  the  extremely  low  density 
mesosphere,  this  filter,  which  strongly  damps  these  modes  while  not 
affecting  the  longer  waves,  is  applied  only  to  the  upper  portion  of 
the  model  during  each  semi -imp! icit  step. 
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3.5  Momentum  Budget 

In  order  to  understand  the  time  evolution  of  zonal  wind,  IT,  we 

need  a  diagnostic  expression  for  3u/3t.  We  start  by  taking  the  zonal 

mean  of  Eq.  (2.7)  of  Holton  (1975)  and  using  the  relationship  between 

basic  density  p  and  scale  height  H  of 
o 

3p 

1  o  =  -1  (192) 

p  W'  T 
o 


to  obtain 


_ _ _ _  2 

_ l] _ [(uv  +  u'v')cos  <)>] 

2  3<J> 

a  cos  <t 


_  3(wiT)  +  wu  _  3(w'u' )  +  wV  +  2  a  V  sin  d>  .  H93) 

__  ^  ___ —  -pp- 


where  u,  v,  and  w  take  on  their  customary  meaning  and  is  the  earth's 
angular  velocity.  By  multiplying  the  above  equation  by  exp[-(z-12)/2H], 
denoted  by  e(z),  we  have 

e(z) OTT/at)  =  FM  +  FB  +  FC  +  FD  +  FE  ,  (194) 


where 


FM  =  e(z)(-  1  9u'v’  +  2u V  tan  »  )  , 

3  <r<j)  3 


(195) 


FB  =  e( z) ( -  1  lii 1  +  2uv  tan  $  ), 
a  3$  a 


(196) 
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FC  =  e(z)2ft 

v  sin 

<t> 

(197) 

FA  =  e(z)(- 

awu  + 

wuT 

“FT 

(198) 

~5z 

FE  =  e  ( z )  (  - 

aw'u' 

+  w'u')  . 

(199) 

az 

nr- 

Hence,  we  can  examine  the  various  contributions  to  the  mean  zonal  wind 
changes  by  horizontal  and  vertical  terms  with  respect  to  both  the  mean 
flow  and  the  eddies. 

For  tropospheric  computations,  we  first  compute  IT,  7,  w,  u'v', 
and  wV  on  sigma  levels  and  then  interpolate  to  find  values  at  1.5, 
4.5,  7.5,  and  10.5  km  in  order  to  mesh  with  the  z  levels  in  the  upper 
portion  of  the  model.  We  then  compute  (194  -  199)  for  these  levels. 

3.6  Heat  Budget 

By  taking  the  time  mean  of  the  thermodynamic  equation  used  by  Lordi 
(1978),  we  can  write 

iX  =  -1  9  C(vT  +  v'T')cos  <j>J 

at  a  cos  4>  a$ 

_  _  2 

-3  (wT  +  w-T1)  -  wHN  +  a(z)(T  -T)  (200) 

az  e 

where  T  is  the  local  departure  of  temperature  from  the  level  mean. 
Multiplying  (200)  by  e(z)  yields 

e(z)(3T/at)  =  VTB  +  VTP  +  WTB  +  WTP  +  HWB  +  HNA  (201) 


-•*> 


where 


VTB  =  -  e( z) 

1 

SvT  cos  $> 

(202) 

a  cos  $ 

3  p 

VTP  =  -  e(z) 

1 

9v'T'  COS  <p 

(203) 

a  cos  <j> 

dip 

WTB  =  -  e(z) 

awT 

az 

(204) 

WTP  =  -  e(z) 

aw'T' 

az 

(205) 

HWB  =  -  e(z) 

_  2 
wHN 

nr 

(206) 

HNA  =  e(z)a(z)(T  -T). 

(207) 

e 


Equation  (201)  gives  us  a  diagnostic  expression  for  aT/at  which  is 

primarily  composed  of  mean  and  eddy  fluxes.  For  the  troposphere,  we 

first  compute  the  fluAes  at  sigma  levels  and  interpolate  to  z  levels 

as  with  the  momentum  budget  equation.  We  then  compute  (201  -  207) 

2 

for  these  levels.  We  also  assume  that  N  given  by  (66)  is  constant 
for  the  troposphere  for  these  diagnostic  purposes  as  is  always  the 
case  with  the  stratosphere. 


CHAPTER  4 


MODEL  RESULTS 

The  cases  were  run  for  a  period  of  40  days  each.  As  mentioned  in 

the  previous  chapter,  the  first  case  included  temperature  profiles  for 

the  stratosphere  and  mesosphere,  which  were  similar  to  those  used  by 

Lordi  (1978).  Case  2  zonal  winds  were  about  50T.  stronger  and  zonal 

temperature  deviations  were  nearly  twice  those  of  the  previous  case. 

The  third  case  was  a  tropospheric  run  using  initial  conditions  in  the 

troposphere  identical  to  those  of  the  other  cases.  For  this  last  case, 

a  lid  was  placed  on  the  troposphere  by  assuming  that  w  =  0  as  an  upper 

R 

boundary  condition,  which  implies  that  the  divergence  contribution  at 
and  below  the  tropopause  from  the  upper  atmosphere  is  also  zero.  The 
same  forcing  was  used  in  each  of  the  three  cases  which  we  will  denote 
Cl,  C2,  and  C3,  respectively.  For  cases  Cl  and  C3,  we  used  a  time  step 
of  30  minutes.  C2  required  a  20  minute  step  because  the  stronger 
flow  caused  limitations  in  meeting  linear  stability  criterion. 

Integrations  were  carried  out  on  the  National  Center  for 
Atmospheric  Research  (NCAR)  CRAY-1  computer.  The  program,  which 
Included  several  output  and  diagnostic  subprograms,  required  approxi¬ 
mately  590  k  of  computer  core  storage.  For  each  10  day  forecast,  Cl, 
C2,  and  C3  required  90,  135,  and  20  minutes  of  computing  time, 
respectively. 

The  time  sections  to  be  presented  in  this  chapter  were  based  on 
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12  hourly  data.  Except  for  polar  projections,  all  tropospheric  repre¬ 
sentations  were  obtained  from  data  interpolated  from  sigma  levels  to 
constant  physical  height  levels.  Polar  plots  for  the  troposphere  depict 
the  values  on  a  constant  sigma  surface.  Latitude-height  sections  of 
geopotential  have  been  density  weighted,  but  time  sections  of  geopoten¬ 
tial  have  not.  Geopotential  amplitude  and  phase  calculations  nave  been 
described  by  Lordi  (1978).  In  this  chapter,  we  will  use  vertical 
cross-sections,  polar  projections,  and  time  sections  to  describe  the 
evolution  of  the  individual  cases,  highlighting  areas  of  similarity  and 
dissimilarity  both  among  the  cases  and  with  observational  data. 

4. 1  Cl  Results 

4.1.1  Evolution  of  zonal  mean  velocity, 
temperature,  and  geopotential 

Latitude-height  secions  of  the  mean  zonal  wind,  mean  zonal  tem¬ 
perature  deviations,  and  mean  zonal  geopotential  deviations  graphically 
illustrate  the  dramatic  changes  that  take  place  during  our  simulated 
sudden  stratospheric  warmings.  For  the  zonal  wind.  Figure  5  depicts 
the  zonal  wind  profile  for  days  10,  20,  30,  and  40.  Comparison  of  day 
10  (Fig.  5a)  with  the  initial  state  (Fig.  2a)  shows  a  slight  weakening 
of  the  polar  night  jet  core  and  slight  displacement  of  the  lower 
extension  of  this  jet  toward  the  pole.  There  is  also  a  minor  decrease 
in  the  amplitude  of  the  tropospheric  jet,  which  should  be  expected 
because  of  the  introduction  of  non-zonal  components. 

Figure  5b  shows  more  dramatic  changes  have  taken  place  by  day  20. 
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The  polar  jet  core  near  65  km  has  weakened  to  about  one-half  of  its 
initial  value  and  a  secondary  core  has  developed  in  polar  latitudes 
around  40  km.  This  feature  has  been  observed  by  Quiroz  et  al .  (1975) 
during  pre-warming  stages  of  actual  sudden  warmings.  This  secondary 
jet  maxima  probably  is  the  atmospheric  response  to  the  squeezing  done 
by  the  northward  moving  easterly  regime  from  the  tropics  to  mid¬ 
latitudes  between  20  and  40  km. 

By  day  30  (Fig.  5c)  the  polar  night  jet  has  weakened  further  as 
the  secondary  core  has  been  completely  eliminated  and  replaced  by 
easterlies  which  extend  from  50  km  to  the  surfarp  at  latitudes  north 
of  70N.  Easterlies  also  extend  from  equator  to  pnlp  ,n  the  lower 
stratosphere  completely  decoupling  the  polar  night  jet  form  the  tropo¬ 
spheric  jet  which  has  intensified  slightly  probably  as  a  response  to 
conserve  angular  momentum  as  suggested  by  Koermer  and  Kao  (1980). 

By  day  40  (Fig.  5d),  the  overall  situation  is  similar  to  day  30. 
However,  the  polar  night  jet  shows  definite  signs  of  recovery  by  inten¬ 
sifying  and  starting  to  push  back  into  the  lower  stratosphere. 
Easterlies  are  also  weakening  both  in  the  polar  stratosphere  and 
troposphere. 

In  order  to  highlight  changes  that  take  place  during  the  reversal 
from  westerlies  to  easterlies,  we  have  also  presented  latitude-height 
sections  for  U  on  days  22,  24,  26,  and  28.  Figure  6a  for  day  22  is 
very  similar  to  Figure  5b,  except  that  at  this  time  the  secondary  jet 
is  as  intense  as  the  upper  core.  By  day  24,  however,  the  secondary  jet 
shows  significant  weakening  with  little  change  elsewhere  as  depicted  in 
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Figure  6b.  This  weakening  trend  continues  on  day  26  (Fig.  6c)  and  the 
easterlies  in  the  lower  stratosphere  at  mid-latitudes  also  weaken. 
Easterlies  have  developed  in  the  polar  troposphere  at  all  levels.  On 
day  23  (Fig.  6d),  the  winds  have  reversed  in  the  lower  stratopshere. 
This  situation  is  very  reminiscent  of  the  easterlies  that  developed 
throughout  such  a  large  vertical  extent  during  the  197b-1977  sudden 
wanning  as  reported  by  Taylor  and  Perry  (1977). 

Although  the  evolution  of  the  zonal  winds  for  this  simulated 
warming  presented  in  Figures  5  and  6  shows  many  of  the  characteristics 
of  both  the  wave  1  and  2  forced  warmings  produced  by  Lordi  (1973), 
there  are  also  several  differences.  One  dissimilarity  is  the  contri¬ 
bution  of  the  easterlies  that  move  from  the  tropics  to  mid-latitudes 
in  the  lower  stratosphere.  In  most  of  Lordi 's  simulations,  these 
easterlies  not  only  push  the  lower  of  the  polar  night  jet  northward 
as  also  happens  in  our  simulations,  but  these  easterlies  also  continue 
to  move  northward  after  the  jet  has  weakened  and  physically  replace 
the  westerlies.  In  our  simulations,  the  easterly  core  after  initially 
displacing  to  mid-latitudes  becomes  quasi-stationary  and  the  easterlies 
in  the  polar  stratosphere  develop  almost  instantaneously  rather  than 
by  advection  into  the  region.  Another  difference  is  the  strength  and 
extent  of  the  easterl ies  in  the  upper  stratosphere  and  mesosphere  that 
resulted  from  Lordi 's  (1978)  cal cul at i ons.  The  pattern  depicted  in 
Figures  5  and  6  show  that  the  mesospheric  westerlies  did  not  reverse 
and  the  strength  of  the  easterlies  that  developed  in  lower  po;a> 
regions  was  about  one-third  of  those  produced  in  s  inui  at  i  •  , 
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(1978).  Possible  reasons  for  these  discrepancies  are  (1)  the  strong 
continuous  forcing  used  by  Lordi  at  the  tropopause  and  (2)  the  forcing 
forcing  by  more  than  one  wave  and  was  non-steady  at  the  boundary,  it 
is  reasonable  to  expect  a  somewhat  weaker  outcome,  especially  at  higher 
levels.  Latitude-height  sections  of  zonal  temperature  deviations 
corresponding  to  the  days  depicted  in  Figures  5  and  6  are  presented  in 
Figures  7  and  8.  Comparing  Figure  7a  with  Figure  3a,  we  can  note  some 
warming  in  the  lower  polar  stratosphere  and  slight  cooling  in  the 
mesospheric  regions  above.  By  day  20,  Figure  7b  shows  this  trend  more 
dramatically.  It  also  shows  stratospheric  cooling  and  mesospheric 
warming  in  equatorial  regions.  The  polar  and  equatorial  regions  of  the 
troposphere  reflect,  to  a  small  degree,  the  temperature  change  charac¬ 
teristics  of  the  stratospheric  regions  above.  Figure  7c  continues  the 
trend  through  day  30  and  more  clearly  depicts  the  cooling  in  the  upper 
stratosphere.  By  day  30,  the  temperature  gradient  in  the  lower  strato¬ 
sphere  has  reversed  with  warmer  temperatures  at  the  pole  than  in  upper 
mid-latitudes,  thermally  supporting  the  easterlies  that  have  developed. 
The  temperature  reversal  continues  through  day  40  (Fig.  7d),  but  is 
confined  to  a  more  shallow  layer  In  the  lower  stratosphere  and  upper 
troposphere.  Cooling  was  beginning  to  take  place  in  the  upper  strato¬ 
sphere  and  mesosphere  in  polar  regions.  Figure  8(a  through  d)  depicts 
the  rapid  temperature  changes  during  the  height  of  the  warming  process, 
when  the  polar  westerlies  are  replaced  by  easterlies. 

During  this  simulated  warming,  polar  stratospheric  temperatures 
rose  by  nearly  408K.  Mesospheric  temperatures  cooled  some  15°K  at  the 
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pole.  Near  the  equator  stratospheric  temperatures  cooled  almost  5°K 
with  the  same  amount  of  warming  in  the  mesosphere.  All  of  these 
features  are  consistent  with  observed  changes  in  actual  warmings  and 
have  been  described  in  Houghton  (1978). 

Figures  9a  through  d  depict  the  density  weighted  latitude-height 
sections  of  geopotential  deviation  of  the  wave  1  component  for  days  10, 
20,  30,  and  40.  From  Figure  9a,  we  can  see  even  by  day  10  the  effect 
of  the  tropospause  in  limiting  the  vertical  propogation  of  wave  energy. 
Energy  does  move  into  the  stratosphere  and  favors  a  more  poleward 
location  than  the  original  tropospheric  forcing  component.  By  day  20 
(Fig.  9b),  wave  1  energy  propogates  even  higher  and  another  surge  of 
geopotential  is  emlnating  from  the  troposphere  into  the  lower  regions 
of  the  polar  stratosphere.  By  day  30  (Fig.  9c),  the  geopotential  wave 
1  amplitude  in  the  upper  region  has  weakened  considerably;  but,  in  the 
lower  regions,  it  has  intensified  and  Is  more  disjointed  from  the 
troposphere.  By  the  final  day  of  the  forecast  period  (Fig.  9d),  the 
wave  amplitude  extends  higher  once  again,  but  the  core  in  the  lower 
stratosphere  has  also  weakened.  The  phase  lines  marked  on  Figures  9a 
and  b  show  increasing  phase  with  height  indicating  warm  air  is  being 
transported  northward  in  the  stratosphere  (Holton,  1975).  This  phase 
relationship  has  altered  significantly  by  day  30  (Fig.  9c).  The 
situation  on  day  40  (Fig.  9d)  also  shows  only  weak  variations  of  phase 
with  height  indicating  nearly  barotropic  conditions. 

Figure  10a  shows  that  on  day  10  wave  2  is  stronger  In  the  lower 
stratosphere  than  wave  1  (Fig.  9a),  but  it  does  not  extend  upward  as 
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far  as  wave  1.  It  also  appears  at  lower  latitudes  and  tends  to  be 
linked  more  closely  with  orographic  forcing.  By  day  20  (Fig.  10b), 
wave  2  Intensifies,  but  changes  little  In  vertical  extent.  There  Is 
also  much  greater  discontinuity  with  the  troposphere.  After  the 
reversal.  Figure  10c  shows  a  much  weaker  wave  2  geopotential  In  the 
stratosphere  which  continues  through  day  40  (Fig.  lOd)  until  wave  2 
has  practically  been  eliminated.  This  follows  from  Charney  and 
Drazin  (1961),  since  the  easterlies  serve  a  critical  layer  preventing 
upward  wave  propagation.  According  to  Houghton  (1978)  when  the  source 
of  wave  energy  Is  cut  off  by  a  critical  level,  the  wave  energy  present 
will  gradually  decay.  Wave  1,  although  affected  to  some  extent  by  the 
critical  level,  probably  has  some  components  with  easterly  phase  speeds 
allowing  it  to  still  propagate  from  the  troposphere  to  the  stratosphere. 
The  phase  relationship  with  height  for  wave  2  closely  parallels  that  of 
wave  1.  Figure  10a  and  b  clearly  indicate  increasing  phase  with  height 
and  warm  air  transport  towards  the  pole.  The  phase  relationship  for 
days  30  and  40  (Figs.  10c  and  d)  Is  not  clearly  defined  and  hence  more 
barotroplc. 

4.1.2  Horizontal  Evolution  of  Geopotential, 

Temperature,  and  Tropospheric  Jets 
Initially,  we  started  with  zonal  mean  fields  for  geopotential, 
temperature  and  velocity.  In  order  to  see  the  meridional  and  longi¬ 
tudinal  development  of  these  fields,  we  will  now  show  the  evolution 

4 

using  polar  projections.  Figure  11a  depicts  geopotential  deviation 
at  40.5  km  after  10  days  of  Integration.  Although  basically  remaining 
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circular,  as  It  was  Initially,  the  vortex  gets  displaces  slightly  from 
the  original  pole  location  and  a  high  develops  downstream  of  the 
larger  orographic  peak.  By  day  20  (Fig.  lib),  a  distinct  wave  2 
pattern  has  developed  with  an  elongated  vortex  still  centered  near  the 
pole.  However,  by  day  30  (Fig.  lie),  the  vortex  has  split  and  a  high 
dominates  the  pole.  Figure  lid  Indicates  the  high  is  weakening  and  a 
less  pronounced  wave  2  pattern  by  day  40.  The  polar  vortex  is  also 
starting  to  regain  its  hold  on  the  polar  regions. 

In  the  troposphere  at  approximately  6.7  km,  the  geopotential 
pattern  becomes  essentially  anchored  by  orographic  forcing  as  depicted 
In  Figures  12a  through  d.  The  vortex  does  elongate  and  split  as  In  the 
stratosphere,  but  the  changes  are  more  subtle  than  In  the  stratosphere. 
We  can  see  that  In  the  troposphere  the  waves  tend  to  be  super-imposed 
on  the  mean  condition  with  greater  wave  activity  during  the  latter 
stages  of  Integration  after  the  warming  process. 

At  40.5  km,  temperature  deviations  respond  similar  to  geopotential 
patterns.  First,  there  Is  a  displacement  of  the  cold  center  from  the 
pole  (Fig.  13a).  Next,  the  displacement  continues  and  temperatures  at 
the  pole  begin  to  warm  as  the  high  Intensifies,  as  shown  in  Figure  13b. 
Thirdly,  the  cold  core  splits  (Fig.  10c)  as  warmer  temperatures  domi¬ 
nate  the  pole.  Finally,  the  warming  subsides  and  conditions  try  to 
return  more  closely  to  a  zonal  mean  state  (Fig.  13d). 

Tropospheric  temperature  deviations  as  shown  in  Figures  14a 
through  c  depict  more  clearly  the  development  of  higher  harmonic  wave 
components  after  evolving  from  a  nearly  zonal  state  on  day  10  to  a 
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more  complex  state  by  day  40.  These  patterns  reflect  to  a  slight 
degree  the  changes  that  have  been  occurring  in  the  lower  stratosphere 


The  tropospheric  jet  first  Intensifies  at  locations  upwind  of  oro¬ 
graphic  forcings  as  depicted  In  Figure  15a.  This  response  is  probably 
due  to  the  increased  convergence  Initially  developing  In  these  regions 
as  topography  Is  gradually  Introduced.  By  day  20  (Fig.  5b),  these  jets 
have  Intensified  proportionally  to  the  height  of  the  topography  and 
have  moved  counterclockwise  until  slightly  downstream  of  the  orographic 
peaks.  Another  more  southerly  positioned  jet  has  developed  upstream  of 
the  weaker  peak.  By  day  30  (Fig.  15c),  the  jet  over  the  higher  peak 
now  remains  essentially  anchored,  but  elsewhere  the  jet  locations  are 
more  transient  and  move  westward,  a  trend  that  continues  on  day  40 
(Fig.  15d) . 

4.1.3  Evolutions  of  Zonal  Mean  Velocity 
Components,  Temperatures,  Geo- 
potential  and  Momentum  and  Heat 
Budgets  at  68.9  N 

Since  we  have  now  given  a  detailed  description  of  Cl  evolution,  we 
want  to  delve  more  closely  Into  the  rapid  changes  that  occur  during  the 
warming  process.  Time-sections  are  an  Ideal  tool  for  studying  these 
changes.  We  will  use  both  height-time  sections  to  study  vertical  varia¬ 
tions  and  latitude-time  sections  to  depict  meridional  variations.  For 
height-time  sections,  68.9°N  was  chosen  as  the  latitude  for  depiction 
since  It  depicts  most  generally  features  at  both  lower  and  higher 
latitudes.  The  7.5  km  level  was  chosen  to  represent  the  troposphere  on 
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latitude-time  sections,  while  40.5  and  67.5  km  were  chosen  to  represent 
the  upper  atmosphere. 

Figure  16  shows  time-height  sections  for  zonal  mean  values  of  u, 
v,  w,  and  T'.  During  the  first  15  days  (Fig.  16a), the  changes  In  u  are 
quite  small  in  contrast  to  the  following  15  day  period  when  the  wind 
reverses.  Figure  16b  shows  two  distinct  southerly  component  peaks 
around  45  to  50  km  on  day  16  and  day  27  followed  by  a  much  weaker 
northerly  peak  on  day  32.  The  two  southerly  peaks  correspond  with  two 
shown  on  Figure  16c.  This  relationship  between  7  and  7  Is  a  very 
Important  Ingredient  of  the  warming  process  and  represents  the  major 
portion  of  an  Indirect  circulation  cell.  As  noted  by  Matsuno  (1971), 
this  cell  consists  of  rising  air  in  polar  regions,  equatorward  trans¬ 
port  at  higher  levels,  and  sinking  motion  at  mid-latitudes.  This  cell 
Is  completed  on  the  lower  end  by  the  eddy  transport  of  sensible  heat 
toward  the  pole.  After  the  zonal  wind  reverses  (Fig.  16a),  Figures  16b 
and  c  Indicate  the  presence  of  a  meridionally  direct  circulation.  The 
temperature  height-time  section  (Fig.  16d)  shows  the  most  rapid  warming 
activity  In  the  stratosphere  during  the  same  periods  when  the  Indirect 
cell  circulations  are  also  most  intense.  Mesospheric  cooling  and 
tropospheric  temperature  changes  are  much  smaller  in  comparison  to 
those  occurring  in  the  stratosphere. 

Figures  17a  through  d  depict  height-time  sections  of  geopotential 
(not  weighted  by  density)  for  zonal  waves  1  through  4,  respectively. 
Wave  1  usually  dominates  the  upper  stratosphere  and  mesosphere  (Fig. 
17a),  whereas  wave  2  (Fig.  17b)  prevails  in  the  lower  stratosphere.  By 


comparing  these  two  wavenumbers,  we  note  that  the  maxlmums  and  mlnlmums 
are  nearly  180°  out  of  phase.  As  mentioned  by  Koermer  and  Kao  (1980), 
this  Indicates  the  non-linear  transfer  of  energy  between  waves  1  and  2. 
Wave  1  shows  a  distinct  oscillation  between  the  peaks  at  60  km  with  a 
period  of  17  days.  Wave  2  has  one  very  strong  peak  occurring  right 
before  the  reversal  and  weakens  considerably  after  day  25.  Tropo¬ 
spheric  peaks  slightly  lead  the  stratospheric  ones.  In  view  of  the 
growth  of  waves  1  and  2  during  the  periods  of  kinetic  energy  decrease 
of  the  zonal  mean  flow  (Fig.  16a),  a  transfer  of  kinetic  energy  from 
the  zonal  mean  flow  to  the  waves  may  occur  after  day  15  of  the 
integration. 

Observational ly  Perry  (1967)  showed  that  wave  3  was  a  source  of 
weakening  for  waves  1  and  2.  By  virtue  of  wave-wave  Interaction,  It 
also  can  help  transfer  energy  between  these  waves.  Our  investigations 
support  these  characteristics  as  shown  In  Fig.  17c.  The  first  peak  for 
wave  3  on  day  21  corresponds  to  a  rapid  weakening  of  wave  1  in  that 
region  and  the  wave  2  Intensification  at  lower  levels  a  few  days  later. 
Conversely,  the  remaining  peaks  correspond  to  wave  2  weakening  and 
precede  wave  1  Intensification  during  the  latter  stages  of  the  warming. 
Wave  4  (Fig.  17d)  remains  relatively  small  and  peaks  when  wave  2 
reaches  Its  maximum.  Hence,  wave  4  plays  a  much  more  minor  role  than 
the  longer  planetary  waves  In  this  simulated  warming.  It  appears, 
therefore,  waves  1,  2,  and  3  form  an  active  triad  in  nonlinear  energy 
transfer.  Geopotential  maxlmums  In  the  upper  atmosphere  tend  to  be 
directly  In  phase  with  maxlmums  in  the  troposphere.  The  amplitude 
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of  the  zonal  wind  also  tends  to  be  Inversely  proportional  to  wave 
amplitude.  These  features  were  noted  in  the  1976-1977  warming  by 
Koermer  and  Kao  (1980). 

At  the  end  of  Chapter  3,  we  introduced  expressions  to  diagnose 
zonal  wind  and  temperature  changes.  By  computing  and  summing  all  the 
terms  on  the  right  side  of  equation  (194),  we  have  the  instantaneous 
rate  of  change  of  zonal  velocity  depicted  as  a  latitude-time  section 
in  Figure  18a.  The  largest  negative  stratospheric  values  Indicate 
decreasing  westerlies  and  occur  on  days  15  and  24.  As  Indicated  by 
16a,  these  westerlies  do  decrease  rapidly  between  these  dates,  making 
these  calculations  consistent.  Instantaneous  velocity  changes  in  the 
troposphere  are  much  more  complex  and  not  as  clearly  defined. 

The  largest  terms  that  make  up  au"/3t  are  FM  and  FC  as  indicated  by 
Figures  18b  and  c.  However,  these  terms  have  opposite  sign  and  nearly 
cancel  open  another  in  the  lower  stratosphere  as  from  day  10  to  day  21. 
During  the  period  from  day  21  to  day  26,  a  period  characterized  by  a 
rapid  weakening  of  the  zonal  wind,  the  FC  term  dominates  in  the  lower 
stratosphere,  resulting  in  westerly  decelerations  from  coriolis  torque, 
which  highlights  the  Importance  of  the  Indirect  cell  circulations  in 
the  warming  process.  During  early  warming  stages,  when  FM  and  FC  are 
nearly  cancelling  one  another,  indicating  the  "non-interaction  effect" 
of  the  meridional  convergence  of  v'u'  on  the  mean  flow,  the  largest 
term  contributing  to  the  decrease  of  the  stratospheric  westerlies  is 
FB,  depicted  in  Figure  18c  with  a  peak  around  day  15  at  40  km.  At 
lower  levels,  the  usually  very  small  FA  term  also  peaks  (Fig.  18d)  and 


also  contributes  to  the  deceleration  of  the  westerlies.  The  FE  term 
(Fg.  18c)  contributes  the  least,  but  also  has  a  negative  peak  around 
day  14  at  30  km  and  does  slightly  aid  in  decreasing  the  westerlies. 

The  momentum  budget  for  Case  1  Indicates  two  distinct  processes 
involved  In  the  deceleration  of  the  polar  night  jet.  During  the  first 
20  days  of  model  Integrations,  an  Indirect  cell  forms,  a  process 
described  by  Matsuno  (1971).  However,  the  convergence  of  westerly 
eddy  momentum  flux  at  higher  latitudes  is  nearly  compensated  by  the 
Coriolis  torque,  basically  negating  the  effects  of  the  indirect  cell 
on  the  rate  of  change  In  the  zonal  flow.  Since  the  larger  terms 
cancel,  the  deceleration  that  occurred  during  this  stage  of  the 
warming  can  be  attributed  to  the  smaller  terms  which  have  their  largest 
negative  values  In  the  stratosphere.  From  day  21  to  26,  the  small 
terms  become  negligible  and  an  imbalance  develops  between  FM  and  FC 
with  FC  slightly  more  dominant.  With  the  negative  term  larger,  rapid 
deceleration  can  take  place,  leading  to  easterlies. 

As  shown  In  Figure  19a,  the  instantaneous  temperature  changes  at 
68. 9N  are  not  as  clearly  defined  in  the  stratosphere  as  were  the  zonal 
wind  changes.  They  are  positive  through  most  of  the  stratosphere, 
indicating  warming  and  negative  in  the  mesosphere  indicating  cooling. 

As  with  the  momentum  budget,  there  are  two  large  terms— YTP  and  HWB— in 
the  zonal  heat  budget  as  indicated  in  Figures  19b  and  c,  respectively. 
However,  they  are  nearly  equal  in  magnitude  and  opposite  in  sign, 
indicating  that  eddy  heat  flux  convergence  is  basically  compensated  by 
adiabatic  cooling  as  the  large  peaks  of  Figures  19b  and  c  are  nearly 
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in  phase  during  the  entire  forecast  period.  This  Indicates  that  there 
exists  a  "non-interaction  effect"  of  the  meridional  convergence  of 
meridional  flux  of  sensible  heat  on  the  vertical  motion.  The  smaller 
terms  depicted  in  Figures  19d  through  g  show  that  among  them  WTB  plays 
the  most  significant  role  In  warming  the  stratosphere,  but  even  it  is 
quite  small.  Perhaps  the  lack  of  definition  of  the  heat  budget  is  due 
to  the  choice  of  latitude  In  this  case,  since  temperatures  change  more 
dramatically  closer  to  the  pole. 

4.1.4  Evolutions  of  Zonal  Mean  Velocities 
Temperature,  Geopotential,  and 
Momentum  and  Heat  Budgets  at  Various 
Heights 

Latitude-time  sections  can  better  show  these  north-south  relation¬ 
ships  and  are  Introduced  here.  Latitude-time  sections  of  zonal 
velocity  at  7.5,  40.5  and  67.5  km  are  presented  in  Figures  20a  through 
c.  Few  changes  occur  in  the  troposphere  except  for  the  nearly  instan¬ 
taneous  switch  to  easterlies  around  day  29  in  polar  regions.  In  the 
stratosphere  at  40.5  km  after  10  days,  we  can  see  the  displacement  of 
the  polar  jet  towards  the  pole  with  a  dramatic  weakening  of  the  jet 
and  ultimate  reversal  between  day  23  and  day  28.  Figure  20b  also 
clearly  shows  that  easterlies  from  mid-latitudes  do  not  displace  pole- 
ward,  but  rather  the  easterlies  develop  on  their  own.  At  high  levels, 
we  can  note  from  Figure  20c  that  higher  latitudes  show  a  southward 
displacement  of  the  polar  night  jet,  the  decreasing  amplitude  of  the 
jet,  and  reintensification  of  the  jet  during  the  later  portion  of  the 
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period  depicted. 

Figures  21a  and  b  clearly  indicate  the  two  major  indirect  cells 
that  develop  during  the  warming  process.  Figure  21a  shows  the  southerly 
component  linking  the  areas  of  strong  vertical  motion  (Fig.  21b). 

Weaker  direct  cells  develop  simultaneously  at  lower  mid-latitudes. 

Two  direct  cells  also  form  in  the  polar  regions  one  with  a  maximum 
amplitude  on  day  22,  when  the  indirect  circulation  at  lower  latitudes 
is  at  a  relative  minimum  and  the  other  during  maximum  easterly 
development. 

Temperature  deviations,  depicted  on  Figures  22a  through  c,  show 
warming  at  7.5,  40.5,  and  67.5  km  in  higher  latitudes  and  cooling  to 
some  extent  at  latitudes  near  the  equator.  Even  though  the  changes  in 
the  troposphere  are  quite  small,  they  reflect  the  changes  that  have 
occurred  in  the  stratosphere. 

Wave  1  geopotential  amplitude  in  the  troposphere  (Fig.  23a)  shows 
a  periodicity  on  the  order  of  10  to  15  days.  The  first  two  peaks  are 
located  around  45°N  where  the  maximum  orographic  forcing  was  positioned. 
The  final  peak  at  the  end  of  the  forecast  period  was  further  north.  At 
40.5  km,  wave  1  amplitude  also  shows  vacillation  with  about  the  same 
periodicity  as  the  troposphere,  but  favors  polar  latitudes  as  indicated 
in  Fig.  23b.  The  first  two  tropospheric  maximums  slightly  lead  the 
stratospheric  maximums,  a  feature  observational ly  consistent  with  data 
presented  by  Labitzke  (1978b)  and  O'Neill  and  Taylor  (1979).  At  even 
higher  levels  (Fig.  23c),  the  wave  amplifies  due  to  decreasing  density 
and  only  two  major  peaks  in  wave  1  amplitude  are  present.  The  first 
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maximum  shows  a  definite  lag  with  lower  level  values,  but  the  second 
one  is  nearly  in  phase  with  that  at  40.5  km. 

Figures  24a  through  c  show  wave  2  geopotential  amplitude  at  the 
same  levels  as  above.  Peaks  at  all  levels  are  generally  not  in  phase 
with  wave  1  peaks  at  the  same  level,  indicating  energy  transfer  between 
these  waves  as  mentioned  by  Koermer  and  Kao  (1980).  Despite  a  vacil¬ 
lating  behavior  in  the  troposphere,  at  40.5  km  the  situation  is  more 
steady  with  essentially  one  major  peak.  At  higher  levels,  the  maximum 
occurs  about  5  days  after  the  40.5  km  maximum,  but  the  values  are  much 
smaller  than  they  were  for  wave  1  at  the  same  level. 

Unlike  waves  1  and  2,  wave  3  (Figs.  25a  through  c)  gets  little  of 
its  energy  from  the  troposphere,  as  we  can  see  by  examining  the  phase 
relations  at  higher  levels  versus  that  at  tropospheric  levels.  Even  in 
the  troposphere,  the  source  of  wave  3  excitation  in  this  model  is  from 
wave-wave  and/or  wave-mean  flow  interactions.  At  40.5  km,  the  wave- 
wave  interaction  is  quite  evident,  since  the  first  peak  corresponds  to 
wave  1  dissipation  and  the  second  peak  to  the  very  large  and  fast  wave 
2  dissipation.  A  similar  relationship  holds  at  67.5  km. 

Wave  4  generally  peaks  with  wave  2  in  the  stratosphere  and 
apparently  leads  wave  4  development  in  the  troposphere,  where  wave  2 
and  4  are  more  out  of  phase  then  in  phase.  As  was  the  case  with  wave 
3,  stratospheric  amplitudes  are  small,  especially  at  high  levels,  as 
we  can  see  from  Figures  26 a  through  c. 

If  we  look  at  the  momentum  budget  at  40.5  km,  we  can  see  more 
clearly  the  physical  processes  involved  in  the  wind  reversal .  Figure 


27a  clearly  shows  the  large  deceleration  taking  place  at  high  latitudes 
from  day  12  through  27.  Again  Figures  27b  and  c  Indicate  that  the  eddy 
momentum  flux  convergence  Is  nearly  balanced  by  the  corlolls  torque 
with  the  latter  being  slightly  larger.  There  Is  little  activity  at 
lower  latitudes.  Of  the  remaining  terms  In  Figures  27d  through  f,  only 
the  mean  momentum  flux  convergence  shows  a  significant  contribution  to 
the  deceleration  of  the  polar  night  jet  and  ultimate  reversal. 

The  heat  budget  at  40.5  km  is  not  clearly  defined,  as  the 
instantaneous  temperature  changes  are  small  and  the  horizontal  eddy 
heat  flux  convergence  is  basically  offset  by  adiabatic  coding,  as 
shown  in  Figures  28a  through  c.  Of  the  four  remaining  small  term 
(Figs.  28d  through  g),  generally  all  but  the  vertical  eddy  heat  flux, 
WTP ,  are  positive  during  the  warming  stages  and  assist  the  warming 
process  at  this  level.  At  lower  levels  (not  shown  here),  horizontal 
eddy  heat  flux,  which  forms  the  bottom  part  of  the  indirect  cell,  is 
more  dominant  and  does  contribute  to  the  warming  of  the  stratosphere. 

4.2  Case  2  Results 

4.2.1  Evolution  of  Zonal  Mean  Velocity, 

Temperature  and  Geopotential 

Latitude-height  sections  of  the  mean  zonal  wind  for  case  C2  with 
stronger  winds  in  the  upper  atmosphere  for  days  10,  20,  30,  and  40  are 
presented  in  Figures  29a  through  d,  respectively.  Up  to  day  20,  the 
evolution  of  these  profiles  are  very  similar  to  Cl.  However,  C2  does 
not  always  maintain  separate  split  cores  as  indicated  in  Figure  29b 
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where  only  the  lower  level  core  remains.  The  tropical  easterlies 
tend  to  displace  the  lower  extension  of  the  polar  night  jet  to  higher 
latitudes,  but  the  critical  line  (u  *  0)  does  not  advance  as  far  or  as 
rapidly  northward  as  in  Cl  (Fig.  5b).  The  easterlies  that  do  develop 
in  the  lower  mid-latitudes  of  the  stratosphere  are  more  intense. 

By  day  30  (Fig.  29c),  a  double  core  is  again  present  and  the 
overall  situation  is  very  similar  to  that  of  day  20  for  Cl  (Fig.  5b). 
The  intensity  of  the  polar  night  at  higher  levels  has  weakened  con¬ 
siderably,  but  in  the  lower  levels,  it  is  still  farily  strong.  Some 
noticeable  differences  between  C2  and  Cl  are  the  downward  propagating 
easterlies  at  subtropical  latitudes  in  the  mesosphere  in  C2  and  the 
more  concentrated  depiction  of  the  tropospheric  jet  in  C2.  By  day  40 
(Fig.  29d),  the  lower  extension  of  the  polar  night  jet  weakened, 
causing  the  secondary  jet  to  vanish,  while  the  jet  at  higher  levels 
in  the  stratosphere  shows  intensification. 

To  highlight  the  most  dramatic  changes  that  occurred  during  C2 
integrations,  we  have  also  presented  the  zonal  wind  profiles  at  2  day 
intervals  between  days  30  and  40.  On  day  32  (Fig.  30a),  the  easterly 
jet  Intensifies  and  pushes  down  into  the  mesosphere.  This  displaces 
the  polar  night  jet  northward  causing  intensification  of  the  jet  in 
the  upper  levels  of  the  stratosphere.  By  day  34  (Fig.  30b),  the 
easterlies  in  the  mesosphere  continued  to  Intensify  as  does  the  polar- 
night  jet  at  higher  latitudes.  By  day  36  (Fig.  30c),  the  mesospheric 
easterlies  have  started  to  subside  and  the  polar  jet  begins  to  return 
to  more  southerly  latitudes.  As  it  regresses,  the  lower  extension  of 


this  jet  weakens  significantly.  This  trend  continues  on  day  38  (Fig. 
30d),  but  apparently  bottoms  out  as  the  jet  In  the  lower  polar  strato¬ 
sphere  shows  signs  of  Intensification  by  day  40  (Fig.  29d).  The 
evolution  of  C2  is  much  more  dramatic  at  higher  levels  than  Cl  and 
indicates  that  strong  winds  serve  as  a  source  to  enhance  wave  activity 
throughout  the  stratosphere  and  mesosphere. 

C2  temperature  deviations  respond  similar  to  Cl,  but  not  as  fast 
as  Cl  through  the  first  20  days  of  Integrations  (Figs.  31a  and  b). 

By  day  30  (Fig.  31c),  there  Is  significant  warming  in  the  mid  and  low 
stratosphere  and  cooling  in  the  upper  stratosphere  and  mesosphere. 

A  slight  temperature  reversal  has  even  occurred  between  40  and  50  km. 
All  of  these  changes  In  polar  regions  are  accompanied  by  equatorial 
cooling  of  the  stratosphere.  By  day  40  (Fig.  31d),  the  temperature 
reversal  is  no  longer  present  in  the  mid-statosphere,  but  lower 
stratospheric  levels  have  warmed  significantly.  Elsewhere,  the 
cooling  trends  are  beginning  to  reverse. 

A  closer  look  between  days  30  and  40,  reveals  that  the  temperature 
reversal  in  the  mid-stratosphere  near  the  pole  was  short  lived  and  not 
even  present  by  day  32  (Fig.  32a).  During  this  ten  day  period  the 
mid-stratosphere  alternately  warmed  and  cooled  as  indicated  by  Figures 
32a  through  d,  but  the  lower  stratosphere  generally  continued  to  warm 
by  nearly  20°K.  However,  by  day  40  (Fig.  31d),  cooling  was  apparent 
at  all  levels  in  the  polar  stratosphere. 

Wave  1  geopotential  (density  weighted)  and  corresponding  phase 
angles  with  height  for  C2  on  days  10  and  20  (Figs.  33a  and  b)  are  quite 


similar  to  those  of  Cl  (Figs.  10a  and  b).  Maximum  amplitudes  for  C2 
are  slightly  south  of  the  locations  In  Cl  and  they  also  extend  to 
slightly  higher  stratospheric  levels.  The  geopotential  extension  from 
the  troposphere  upward  and  northward  Into  the  polar  stratosphere  Is 
similar  for  both  cases,  but  the  amplitude  of  C2  Is  larger.  By  day  30 
(Fig.  33c),  wave  1  develops  a  tremendous  amplitude  thorughout  the 
entire  polar  stratosphere,  a  feature  not  present  In  Cl.  This  large 
peak  is  nearly  disjointed  form  the  troposphere.  Figure  33  shows  a 
similar  pattern  existing  again  on  day  40,  indicating  a  strong  wave  1 
component.  As  with  Cl,  phase  angles  for  C2  are  more  baroclinic  at 
early  stages  of  the  warming  and  become  essentially  barotropic  at  later 
stages. 

Figures  34a  through  d  show  wave  2  geopotential  evolution  for  C2. 
The  phase  relationship  evolves  similar  to  that  mentioned  above  for  wave 
1  with  the  greatest  increase  of  phase  with  height  on  day  10.  Wave  2 
also  penetrates  higher  levels  when  the  stronger  westerlies  are  present 
with  amplitudes  slightly  south  of  the  Cl  locations  in  the  stratosphere 
and  with  C2  amplitudes  spread  over  a  broader  area.  Wave  2  weakening 
In  the  stratosphere  during  the  latter  phases  of  the  integration  period 
are  quite  similar  for  both  C2  and  Cl,  despite  the  fact  that  the  winds 
did  not  reverse  In  C2.  This  may  Indicate  the  presence  of  a  strato¬ 
spheric  vacillation  cycle  as  discussed  by  Holton  and  Mass  (1976). 
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4.2.2  Horizontal  evolution  of 

Geopotential,  Temperatures, 
and  Tropospheric  Jets 

Polar  projections  of  geopotential  deviatons  at  40.5  km  show 
little  disfiguration  of  the  polar  vortex  of  day  10  (Fig.  35a).  By  day 
20  (Fig.  35b),  high  pressure  cells  do  develop  and  tend  to  squeeze  the 
vortex  In  a  wave  2  pattern  similar  to  Cl  (Fig.  11b).  However,  the 
vortex  continues  to  hold  in  polar  regions  on  day  30  (Fig.  35c)  with 
high  pressure  cells  no  longer  on  opposite  sides  of  the  vortex,  but 
tending  to  merge,  which  they  effectively  do  by  day  40  (Fig.  35d). 

Acting  together,  the  polar  vortex  gets  displaced,  but  not  far  enough 
for  easterlies  to  dominate  the  polar  regions. 

In  the  troposphere,  the  evolution  of  geopotential  for  C2  during 
the  first  half  of  the  integration  period  (not  shown  here)  is  nearly 
identical  to  that  for  Cl  (Figs.  12a  and  b).  During  the  second  half 
(Figs.  36a  and  b),  C2  does  show  some  subtle  differences  with  Cl  (Figs. 
12c  and  d),  with  C2  having  a  more  pronounced  wave  3  and  4  pattern  than 
Cl.  Usually,  the  geopotential  configuration  near  the  tropopause  at 
13.5  km  was  very  close  to  that  of  7.5  km.  However,  during  the  last  ten 
day  period  of  C2  when  warming  activity  was  the  strongest,  a  strong 
wave  3  pattern  developed  (Figs.  36c  and  d),  which  is  often  associated 
with  blocking  in  the  troposphere.  This  feature  Is  not  as  promlnant  at 
lower  levels  because  of  the  dominance  of  waves  1  and  2  from  orographic 
forcing. 

The  C2  plots  of  temperature  deviation  for  40.5  km  on  days  30  and 
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40  and. the  corresponding  7.5  km  results  are  shown  on  Figures  37a 
through  d.  Patterns  during  earlier  stages  of  the  Integrations  were 
similar  to  those  of  Cl  (Figs.  13a  and  b).  In  comparison  with  Cl  at 
later  stages  (Figs.  13c  and  d),  we  can  see  that  wave  1  Is  much  more 
promlnant  In  the  stratosphere  for  C2.  The  rapid  migration  of  the  cold 
core  about  the  pole  for  this  case  Indicates  large  traveling  components 
and  wave  transience.  In  the  troposphere,  patterns  are  still  most 
strongly  Influenced  by  orography  and  hence  are  quasl-statlonary,  but 
C2  (Figs.  37b  and  d)  shows  more  diversity  in  wave  activity  than  Cl 
(Figs.  14b  and  d). 

The  tropospheric  jet  maxima  for  Cl  (Figs.  15a  through  c)  were 
quite  similar  to  C2  (not  shown).  However,  the  C2  maxima  were  generally 
weaker,  tending  to  make  the  overall  flow  more  zonal  than  in  Cl.  The 
apparent  anchoring  of  one  jet  to  the  highest  topography  was  present  in 
both  cases,  as  was  the  transience  of  the  maxima  elsewhere. 

4.2.3  Evolutions  of  Zonal  Mean 

Velocity  Components,  Temperature, 

Geopotential  and  Momentum  and 
Heat  Budgets  at  68.9  N 

Even  though  a  major  warming  according  to  WMO  standards  failed 
to  develop  for  C2  as  it  did  for  Cl,  the  so-called  "minor"  warming  was 
even  more  dynamic.  The  time  section  of  zonal  wind  (Fig,  38a)  shows 
strong  vacillations  after  day  20  with  a  periodicity  of  around  7  days. 
During  the  pre-warming  stages,  the  jet  actually  Intensifies.  The 
zonal  mean  v  (Fig.  38b)  and  w  (Fig.  38c)  also  Indicate  5  Indirect  cell 
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circulations  with  the  most  Intense  cell  around  day  25  corresponding 
to  the  period  of  greatest  wanning  as  shown  In  Figure  38d.  The  period¬ 
icity  of  the  first  two  of  these  cells  Is  quite  similar  to  that  of  Cl, 
but  the  C2  circulations  are  more  Intense,  especially  the  second  one. 

The  smaller  Indirect  circulations  that  developed  after  this  large  one 
were  unique  to  C2,  since  In  Cl  a  direct  circulation  developed  after 
the  easterlies  developed  In  the  polar  stratosphere.  At  68.9°N,  mid- 
stratospheric  temperatures  peaked  around  day  29.  It  was  preceded  by 
rapid  warming  and  followed  by  fairly  rapid  cooling  as  indicated  In 
Figure  38d.  Warming  of  the  lower  stratosphere  took  place  gradually 
almost  from  the  start  of  integrations,  but  was  most  Intense  from  days 
31  to  35. 

As  shown  in  Figure  39a,  wave  1  geopotential  for  C2  clearly 
indicates  stronger  and  more  rapid  vacillations  than  occur  In  Cl  (Fig. 
17a)  during  the  second  half  of  the  Integration  period.  During  these 
stages  when  C2  Is  most  active,  wave  1  peaks  are  exactly  In  phase  with 
if  minimums  (Fig.  38a),  indicating  that  a  portion  of  wave  1  energy  is 
being  extracted  form  the  mean  flow.  This  Is  similar  to  the  observed 
data  presented  by  Koermer  and  Kao  (1980).  Because  of  the  strong 
fluctuations  of  wave  1  (Fig.  39a)  and  from  the  polar  representation 
of  geopotential  (Figs.  35a  and  b).  It  Is  apparent  that  the  polar  vortex 
Is  wobbling  like  a  spinning  top  about  the  pole.  When  it  becomes 
centered  over  the  pole,  the  amplitude  of  wave  1  suddenly  decreases. 

When  It  shifts  away  from  the  pole,  there  Is  rapid  intensification  of 
wave  1.  Generally  tropospheric  peaks  slightly  lead  stratospheric 
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maxima. 

Wave  2  geopotential  {Fig.  39b)  for  C2  Is  also  stronger  and 
varies  more  frequently  during  the  latter  stages  of  the  forecast 
period  than  It  was  for  Cl,  reaffirming  the  Importance  of  the  critical 
level  that  develops  when  the  easterlies  extend  from  equator  to  pole  In 
cutting  off  the  source  of  wave  excitation  In  the  stratosphere.  Wave  2 
Is  essentially  out  of  phase  with  wave  1,  Indicating  (as  In  Cl)  the  non¬ 
linear  transfer  of  energy.  This  transfer  is  again  aided  by  wave  3 
(Fig.  39c),  which  intensifies  more  slowly  at  first  than  it  did  In  Cl; 
but  then  Intensifies  and  oscillates  rapidly  during  the  latter  stages  of 
the  forecast  period.  As  with  Cl,  wave  3  maxima  for  C2  tend  to  occur 
during  the  periods  when  wave  1  is  Increasing  and  wave  2  is  decreasing 
or  vice-versa.  Wave  4  (Fig.  39d)  is  essentially  In  phase  with  wave  2 
but  is  of  much  higher  amplitude  than  it  was  for  Cl  (Fig.  17d). 

The  strong  stratospheric  vacillations  indicated  by  geopotential 
oscillations  of  the  wave  components  in  C2  do  not  show  In  the  analysis 
of  Holton  and  Dunkerton  (1978).  Using  a  linearized  model,  they  suggest 
that  waves  which  are  able  to  propagate  to  large  heights  produce  weaker 
amplitude  vacillations  than  those  that  can  propagate  only  to  lower 
levels.  Cl  corresponds  to  this  latter  case  In  that  waves  do  not 
propagate  to  high  levels  when  the  zonal  winds  are  weak.  C2  with 
stronger  winds  and  consequently  higher  wave  propagation  corresponds 
to  the  former.  However,  C2  produced  the  largest  vacillations  In  our 
Integrations  compared  to  Cl.  This  difference  is  probably  due  to  the 
non-linear  effects  and  due  to  the  presence  of  four  waves  In  our  model. 


The  strong  zonal  winds  of  C2  provide  a  large  energy  source  from  which 
all  waves  can  extract  energy.  If  the  waves  are  stronger,  the  wave-wave 
Interaction  transfers  of  energy  will  probably  be  larger  and  greater 
vacillations  result  along  with  greater  wave  transience.  However,  If 
only  the  waves  and  mean  flow  Interact,  this  may  not  be  the  case,  as 
shown  by  Holton  and  Dunkerton  (1978).  This  discrepancy  clearly  shows 
the  Importance  of  including  wave-wave  interactions  in  the  modeling  such 
phenomena  as  stratospheric  warmings. 

Instantaneous  velocity  changes  In  the  stratosphere  are  also  more 
pronounced  for  C2  as  Indicated  by  Figure  40a.  As  should  be  expected, 
this  means  that  the  momentum  budget  terms  should  be  more  pronounced  to 
cause  the  larger  changes,  primarily  In  decelerating  the  zonal  wester¬ 
lies  of  the  stratosphere.  As  with  Cl,  the  large  eddy  momentum  flux 
convergence  (FM)  (Fig.  40b)  Is  compensated  to  a  large  extent  by  the 
cor  t1 1 s  torque  (FC)  (Fig.  40c).  During  the  early  stages  of  C2  inte¬ 
grations,  FM  (which  peaks  on  day  15  of  this  period)  is  slightly  larger 
than  FC,  aiding  the  acceleration  of  the  westerlies  as  Indicated  on 
Figure  38a.  During  the  next  large  oscillation,  the  peaks  nearly 
coincide  and  cancel  because  of  their  opposite  sign.  During  the  latter 
stages  of  the  forecast  period  marked  by  frequent  vacillation,  these 
terms  change  rapidly  and  do  not  always  balance  one  another,  as  Is  the 
case  between  days  35  and  37,  where  FC  Is  more  dominant.  This  Imbalance 
probably  Is  most  responsible  for  the  final  deceleration  phase  of  the 
westerlies  as  shown  In  Fig.  38a.  The  smaller  terms  (Figs.  40d  through 
f)  show  two  major  peaks  In  the  stratosphere  corresponding  to  the  first 
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two  maxima  of  FM  and  FC  (Figs.  40b  and  c)  and  closely  linked  to  the 
Indirect  cell  mechanism.  FB,  FA,  and  FE  all  tend  to  decelerate  the 
stratosphere  westerlies  at  higher  latitudes  when  they  are  large  in 
amplitude.  Since  they  are  still  often  small  In  relation  to  the  large 
terms  (FM  and  FC),  it  makes  their  contribution  more  dubious.  For 
example,  the  zonal  westerlies  in  the  stratosphere  still  increased 
during  the  pre-warming  stages  despite  the  fact  that  vertical  fluxes  and 
mean  momentum  flux  were  working  against  this  since  there  was  sufficient 
imbalance  between  FM  and  FC  to  result  in  acceleration.  However,  these 
terms  can  play  a  role,  as  during  the  second  major  cycle  when  FM  an  FC 
effectively  cancel  and  zonal  wind  changes  still  occur.  FB,  FA,  and  FE 
contribute  little  during  the  latter  stages  of  the  forecasts. 

The  net  result  of  instantaneous  temperature  changes  (Fig.  41a) 
is  not  clearly  defined  because  of  the  cancelling  effect  of  the  contri¬ 
buting  terms.  From  Figure  41b,  we  can  see  the  lower  portion  of  the 
indirect  cells  alluded  to  earlier.  This  eddy  meridional  heat  transport 
is  counteracted  to  a  large  extent  by  the  adiabatic  cooling  effects  of 
HWB  (Fig.  41c).  During  the  first  pulse,  the  net  heating  between  days 
10  and  18  in  the  stratosphere  as  indicated  on  Figure  38d  may  be  due  to 
WTB  (Fig.  41e)  since  VTP  and  HWB  nearly  cancel.  However,  WTP  (Fig. 

41 f)  reduces  the  extent  WTB  and  the  other  small  terms  VTB  (Figs.  41 d 
and  f)  show  little  contribution.  During  the  second  warming  pulse 
between  days  24  and  29  (Fig.  38d),  YTP  (Fig.  41b)  with  weak  assistance 
from  WTB  (Fig.  41e)  apparently  dominates  HWB  (Fig.  41c)  and  the  other 
small  terms,  resulting  in  mid-stratospheric  warming.  From  days  20  to 
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33,  this  trend  reverses  and  cooling  takes  place  (Fig.  38d). 

4.2.4  Evolutions  of  Zonal  Mean 
Velocities,  Temperature, 

Geopotential,  and  Momentum  and 
Heat  Budgets  at  Various  Heights 

Latitude-time  sections  for  zonal  wind  are  shown  in  Figures  42a 
through  c.  Tropospheric  wind  evolution  (Fig.  42a)  is  similar  to  that 
of  Cl  (Fig.  20a),  except  easterlies  did  not  develop  for  C2.  At  40.5  km 
(Fig.  42b),  the  slower  evolution  of  C2  is  evident  although  the  poleward 
displacement  of  the  jet  is  similar  to  Cl  (Fig.  20b).  The  time  evolution 
of  mean  zonal  wind  at  67.5  km  (Fig.  42c),  the  situation  is  quite  dif¬ 
ferent  from  Cl  (Fig.  20c)  with  C2  showing  greater  oscillations  along 
with  more  rapid  decelerations  at  lower  latitudes.  This  is  probably 
the  result  of  easterly  jet  development  in  the  mesosphere  above  these 
latitudes. 

The  major  indirect  cell  circulations  are  clearly  evident  at  higher 
latitudes  from  Figure  43a  and  b.  Both  cells  develop  at  about  the  same 
time  as  they  had  done  during  C2  integrations,  probably  indicating  an 
inherent  vacillation  process  as  proposed  by  Holton  and  Mass  (1976) 
since  at  this  time  waves  were  primarily  extracting  energy  from  the  mean 
flow.  Thermally  direct  circulations  also  simultaneously  develop  in 
mid-latitudes  as  they  had  done  in  Cl.  All  of  these  circulation  cells, 
although  similar,  lie  further  south  for  C2  than  In  Cl. 

Tropospheric  temperature  changes  (Fig.  44a)  for  C2  are  similar  to 
Cl  (Fig.  22a)  except  for  the  sudden  cooling  in  the  polar  regions  that 
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occurred  after  the  initial  warming  during  C2  integrations.  At  40.5  km 
(Fig.  44b),  sudden  warming  in  polar  regions  occurs  from  days  24  to  30, 
followed  by  cooling  almost  as  large  from  days  30  to  34.  At  lower 
latitudes,  the  evolution  is  remarkably  similar  to  Cl  (Fig.  22b).  The 
warming  and  cooling  surges  at  67.5  km  (Fig.  44c)  basically  mimic  those 
of  40.5  km  just  described.  C2  temperature  changes  were  slower  in 
coming  about  than  in  Cl,  but  when  they  did  occur  the  changes  were  much 
more  rapid  and  about  the  same  order  of  magnitude  as  in  Cl. 

In  the  troposphere,  C2  wave  1  geopotential  (Fig.  45a)  looks  much 
like  that  of  Cl  (Fig.  23a)  for  the  first  30  days.  However,  after  the 
wind  reverses  in  Cl,  wave  1  tropospheric  geopotential  appears  much 
larger  in  the  polar  regions.  At  40.5  km  (Fig.  45b),  C2  wave  1 
geopotential  is  essentially  similar  to  that  of  Cl  (Fig.  23b)  with  a 
maximum  amplitude  during  these  pre-warming  stages  slightly  larger  and 
about  10°  further  south  than  Cl.  The  last  13  days  of  each  case  are 
considerably  different  with  tremendous  wave  1  activity  for  C2  as 
described  earlier  in  this  section.  At  even  higher  levels  (Fig.  45c), 
wave  1  appears  to  be  out  of  phase  and  lagging  the  40.5  km  level. 

Higher  periodicities  are  characteristic  of  the  latter  part  of  C2  inte¬ 
grations.  It  Is  during  the  times  of  sudden  wave  1  amplification  that 
westerlies  tend  to  decelerate  most  rapidly. 

Wave  2  geopotential  for  C2  (Fig.  46a)  at  7.5  km  has  the  same  basic 
features  as  Cl  (Fig.  24a),  showing  only  slightly  larger  amplitudes  at 
higher  latitudes.  At  40.5  km  (Fig.  46b),  wave  2  is  out  of  phase  with 
wave  1  (Fig.  45b).  This  is  especially  vivid  from  day  27  to  day  31, 
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where  wave  2  decreased  rapidly  as  wave  1  simultaneously  amplified. 
Recall  that  the  warming  essentially  peaked  during  this  period.  Wave  2 
(Fig.  46c)  is  much  stronger  and  more  active  at  67.5  km  than  it  was  for 
Cl  (Fig.  24c).  Again,  it  shows  an  out  of  phase  relatonship  with  wave  1 
(Fig.  45c). 

Although  wave  3  (Fig.  47a)  evolves  in  the  troposphere  at  about  the 
same  time  as  it  does  in  Cl  integrations  (Fig.  25a),  we  can  see  that  the 
C2  development  is  much  more  intense,  reflecting  the  wave  3  pattern 
discussed  earlier  that  occurs  during  the  latter  stages  of  the  forecast 
period.  At  higher  levels  (Figs.  47b  and  c),  wave  3  is  also  stronger 
than  in  the  previous  case  and  its  vacillations  correspond  to  wave  1  and 
2  energy  transfers. 

Wave  4  geopotential  for  C2  (Fig.  48a)  has  about  the  same  tropo¬ 
spheric  amplitude  as  it  did  in  Cl  (Fig.  26a).  Stratospheric  values 
(Figs.  48b  and  c)  are  in  contrast  much  larger  for  C2  indicating  the 
greater  energy  source  available  initially.  As  with  Cl,  wave  4  maximums 
are  closely  linked  to  wave  2  peaks. 

Latitude-time  sections  for  the  momentum  budget  at  40.5  km  are  pre¬ 
sented  in  Figures  49a  through  f  and  essentially  confirm  the  processes 
described  earlier  based  on  time-height  sections.  FM  is  basically  the 
only  term  trying  to  maintain  or  accelerate  the  westerlies.  FC  provides 
the  greatest  source  for  deceleration  and  in  fact  the  greatest  deceler¬ 
ation  usually  occurs  when  FC  Is  clearly  dominant  over  FM.  During  the 
last  12  days,  the  C2  warming  remains  characterized  by  oscillations 
reflected  in  the  momentum  budget  terms  and  which  is  probably  a  result 
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of  the  rapid  energy  transfer  between  waves. 

The  instantaneous  temperature  change  (Fig.  50a)  also  indicates  the 
oscillating  behavior  of  C2,  especially  at  high  latitudes.  Most  of  the 
fluctuations  can  be  attributed  to  small  imbalances  between  VTP  (Fig. 

50b)  and  HWB  (Fig.  50c).  Warming  is  aided  to  a  small  extent  at  high 
latitudes  by  VTB  and  WTB,  indicated  in  Figures  50d  and  e.  At  somewhat 
lower  latitudes,  WTP  (Fig.  50 f )  contributes  to  cooling.  HNA  contri¬ 
butes  little  to  the  heating  and  cooling  processes  at  this  level. 

In  summary,  C2  is  similar  to  Cl  during  the  first  half  of  the 
forecast  period,  but  C2  becomes  more  complex  during  the  latter  stages 
when  it  becomes  highly  oscillatory  in  all  respects.  During  these  sudden 
changes,  imbalances  in  indirect  circulations  develop  resulting  in 
deceleration  or  acceleration  of  the  zonal  wind  and  heating  or  cooling. 
These  processes  are  in  good  agreement  with  the  basic  mechanism  proposed 
by  Matsuno  (1971).  The  oscillations  themselves  tend  to  be  the  result 
of  wave  transience  and  wave-wave  interactions.  In  this  case  as  in  Cl, 
large  terms  in  the  momentum  and  heat  budgets  nearly  cancel.  However, 
the  slight  imbalances  between  these  large  terms  can  control  the  warming 
process.  It  is  only  when  these  large  terms  essentially  cancel,  that 
the  smaller  terms  become  important  enough  to  affect  significant  wind 
and  temperature  changes.  Even  though  the  warming  generated  in  C2  over 
the  40  day  integration  period  cannot  be  classified  as  major,  it  was  as 
dynamic  as  the  major  Cl  warming  in  its  evolution.  C2  results  also  tend 
to  confirm  the  hypothesis  of  Koermer  and  Kao  (1980)  that  a  strong  polar 
vortex  can  inhibit  major  stratospheric  warmings. 


4.3  Case  3  Results 

Case  3  (C3)  involved  running  the  model  for  the  troposphere  only 


using  the  same  Initial  conditions  for  these  levels  as  were  used  in  Cl 
and  C2.  Identical  forcing  was  also  used.  A  lid  was  placed  on  the 
model  at  12  km  negating  interactions  with  the  stratosphere.  By  doing 
this,  we  are  able  to  highlight  interaction  effects  form  the  strato¬ 
sphere  to  the  troposphere  of  Cl  and  C2  by  comparing  their  tropospheric 
results  with  those  of  C3. 

Latitude  height  sections  of  zonal  velocity  for  C3  are  shown  in 
Figures  51a  through  e.  Relatively  few  changes  from  the  original  pro¬ 
file  (Fig.  51a)  occur  until  late  in  the  integration  period.  At  this 
time  (Figs.  51 f  and  e)  a  dual  jet  develops  in  the  subtropics  and  at 
mid-latitudes.  This  feature  did  not  develop  in  Cl  and  C2,  but  may  be 
the  result  of  upper  boundary  reflection. 

Polar  plots  of  geopotential  deviation  at  approximately  6.7  km  show 
a  very  monotonous  pattern  in  Figures  52a  through  d.  The  wave  activity 
other  than  that  being  forced  by  orography  is  much  less  pronounced  than 
that  which  resulted  by  Inclusion  of  an  upper  atmosphere.  The  tempera¬ 
ture  deviations  (Figs.  53a  through  d)  also  follow  this  trend  with  a 
zonal  type  pattern  throughout  the  forecast  period  in  contrast  with  the 
non-zonal  nature  of  the  temperature  distibutlon  of  Cl  (Figs.  14c  and 
d)  and  of  C2  (Figs.  37c  and  d)  recorded  at  later  stages  of  the 
Integrations.  The  C3  zonal  wind  maxima  (Figs.  54a  through  d)  shows  a 
rather  slow  development  versus  the  other  cases.  There  is  also  no 
apparent  anchoring  of  a  jet  to  topography  in  C3  as  was  noticed  in  Cl 
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and  C2.  Instead,  all  the  maxima  move  around  a  latitude  circle  with  the 
mean  flow. 

At  7.5  km  (Fig.  55a),  changes  in  the  C3  zonal  wind  are  slower, 
weaker  and  steadier  than  the  other  cases.  Polar  easterlies  do  develop 
as  they  did  in  Cl,  but  the  development  is  more  gradual  than  the  sudden 
change  noted  in  Cl  (Fig.  20a).  Temperature  (Fig.  55b)  does  not  signi¬ 
ficantly  change  during  the  entire  forecast  period  in  stark  contrast  to 
the  development  of  Cl  (Fig.  22a)  and  C2  (Fig.  44a).  This  difference 
probably  indicates  an  important  interaction  between  the  stratosphere 
and  the  troposphere  in  that  the  troposphere  reflects,  to  a  small 
degree,  the  stratospheric  temperature  changes.  As  noted  by  Namais 
(1978),  the  mean  temperatures  at  higher  latitudes  in  January  1977  were 
much  warmer  than  normal  following  a  major  sudden  warming  very  similar 
to  that  produced  in  Cl.  Although  many  areas  were  warmer,  some  areas  of 
the  United  States  experienced  the  coldest  January  ever  known  (Wagner, 
1977).  Hence,  the  stratosphere  may  influence  the  troposphere  signi¬ 
ficantly  enough  to  produce  long-term  climatic  anomalies. 

Finally,  the  geopotential  wave  amplitudes  for  waves  1  through  4 
(Figs.  56a  through  d)  show  much  less  variability  for  C3  than  was  shown 
by  the  other  cases.  Although  the  amplitudes  of  wave  1  an  2  are  similar 
and  evolve  similar  to  that  of  Cl  and  C2,  waves  3  and  4  simply  do  not 
develop  with  anywhere  near  the  intensity  as  shown  in  the  previous  cases. 
The  former  item  we  would  expect,  since  waves  1  and  2  in  the  troposphere 
are  controlled  primarily  by  the  orographic  forcing,  which  is  the  same 
for  both  cases.  The  latter  item  may  be  one  reason  why  many  general 


circulation  models  have  their  greatest  problems  In  accurately 
predicting  the  evolution  of  the  planetary  waves  (Sommervllle,  1979) 


CHAPTER  5 


CONCLUSION 

In  this  thesis,  we  have  formulated  a  new  model  to  examine  the  pro¬ 
cesses  Involved  in  a  sudden  stratospheric  warming.  We  have  simulated 
a  major  and  a  minor  warming  during  forty  day  time  integrations.  We 
have  looked  at  the  momentum  and  heat  budget  contributions  and  the  wave 
evolutions  associated  with  these  warmings.  Finally,  we  have  shown 
results  indicating  the  importance  of  including  the  upper  atmosphere  in 
modeling  efforts. 

The  model  presented  is  Innovative  and  can  serve  as  a  building 
block  for  even  more  sophisticated  studies.  Most  spectral  models  that 
include  a  troposphere  in  sigma-coordinates  have  been  based  on  pr  =  0 
(e.g.  Bourke,  1974;  Hoskins  and  Simmons,  1975).  Although  this 
simplifies  the  resulting  equations,  it  does  not  allow  for  enough 
definition  in  the  upper  atmosphere.  By  using  two  separate  coordinate 
systems  and  pr  t  0,  we  can  have  detailed  resolution  of  the  stratosphere 
and  mesosphere  using  log-pressure  coordinates  and  still  use  the  sigma- 
coordinates  for  the  troposphere.  These  latter  coordinates  have 
excellent  terrain  handling  properties  (Kasahara,  1974)  and  they  also 
simplify  the  specification  of  bottom  boundary  conditions. 

The  spectral  form  of  the  model  lets  us  use  semi-implicit  time 
differencing,  which  in  turn  allows  us  to  use  larger  time  steps. 
Spherical  harmonics  also  allow  us  to  perform  horizontal 
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differentiations,  Integrations,  and  fast  Fourier  transforms  exactly  up 
to  the  point  of  truncation.  The  model  Is  also  designed  to  conserve 
mass  and  energy  by  using  appropriate  vertical  differencing  schemes. 

This  and  the  semi-lmpl icit  time  differencing  give  the  model  the 
capability  to  perform  long  term  Integrations. 

Currently,  the  biggest  drawbacks  of  this  model  are  the  computer 
core  size  required  and  the  large  amount  of  computer  time  required  to 
perform  the  time  integrations.  There  are  several  areas  that  could  be 
trimmed  to  reduce  both  of  these  requirements.  The  most  important  would 
be  in  divorcing  various  diagnostic  and  output  routines  from  the  main 
program.  They  could  be  run  separately,  independent  of  the  main  program. 
Another  area  that  could  be  examined  to  promote  greater  efficiency  is 
the  requirement  for  26  levels  In  the  upper  atmosphere.  Experimentation 
might  show  that  several  levels  could  be  eliminated,  which  would  reduce 
both  core  size  and  save  computational  time.  Greater  efficiency  might 
also  be  achieved  by  program  consolidation  which  would  allow  for  array 
reduction. 

The  model  has  great  potential  for  future  studies.  If  greater 
efficiency  can  be  achieved,  the  model  can  be  expanded  in  several  ways. 
The  most  obvious  expansion  would  be  the  Inclusion  of  higher  harmonics 
In  the  troposphere.  Currently,  the  model  is  truncated  at  wave  4  and 
does  not  allow  for  the  Inclusion  of  the  very  Important  synoptic  scale 
waves.  The  additional  harmonics  would  let  energy  cascade  to  these 
meteorologically  important  wavenumbers.  The  stratosphere  Is  probably 
sufficiently  covered,  as  we  can  see  by  small  wave  4  amplitudes  that 
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developed  In  cases  Cl  and  C2,  and  would  not  have  to  be  expanded 
further.  Another  area  of  possible  expansion  would  involve  not  assuming 
symmetry  across  the  equator  and  consequently  going  to  a  true  global 
model.  This  would  allow  for  inclusion  of  topography  representative  of 
both  hemispheres  and  also  Include  transports  across  the  equator.  A 
third  major  area  of  future  research  would  be  the  inclusion  of  more 
physics,  such  as  moisture  and  differential  heating  parameterization. 

All  of  these  modifications  would  add  realism  and  subsequently  even 
greater  accuracy  to  the  model  results. 

In  our  Integrations  of  Cl  and  C2,  we  have  demonstrated  how  sudden 
warmings  can  be  initiated  by  orographic  forcing.  We  have  also  shown 
how  the  evolution  of  these  warmings  Is  dependent  upon  the  Initial 
intensity  of  the  polar  night  jet.  Given  the  same  degree  of  forcing, 
the  weaker  jet  yields  more  easily  to  the  forcing  from  below  and  this 
situation  can  more  readily  lead  to  a  major  sudden  stratospheric  warming. 
A  stronger  jet  is  at  first  less  affected  by  forcing,  but  as  the  waves 
begin  to  extract  energy  from  the  large  amount  available  in  the  mean 
flow,  the  warming  can  become  very  active.  Until  this  point  was  reached, 
however,  the  Initial  evolution  of  Cl  and  C2  was  quite  similar.  Cl 
later  showed  smaller  amplitude  waves  with  longer  periodicity,  especially 
at  higher  levels  versus  the  larger  and  higher  frequency  characteristics 
of  C2. 

The  stratospheric  vacillations  appear  to  be  caused  by  the  process 
of  energy  transfer  resulting  from  wave  to  wave  interactions.  Although 
wave  an  mean  flow  interactions  may  also  account  for  a  certain  degree  of 
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vacillation,  especially  during  the  early  stages  of  our  integrations.  It 
cannot  account  for  the  much  higher  frequency  oscillations  produced  in 
C2.  By  comparing  wave  growth  and  decay  among  the  planetary  waves,  the 
importance  of  the  wave  to  wave  Interactions  to  the  vacillalon  process  is 
quite  evident.  Because  of  the  inability  to  sustain  a  process  for  a  long 
enough  period  in  a  highly  vacillating  regime  with  significant  wave 
transience,  the  C2  warming  most  probably  did  not  become  major. 

Our  integrations  certainly  support  the  indirect  cell  mechanism  as 
outlined  by  Matsuno  (1971)  as  the  prime  Ingredient  In  the  warming 
process.  Although  the  large  terms  Involved  In  momentum  and  heat 
budgets  tend  to  cancel,  a  slight  Imbalance  between  them  can  cause 
large  changes  in  zonal  wind  and  temperature.  Sometimes  they  do  receive 
support  from  the  smaller  contributors  to  these  budgets,  which  become 
increasingly  important  when  the  larger  terms  effectively  negate  one 
another. 

It  is  still  difficult  to  conclude  whether  wave  1  as  supported  by 
Labitzke  (1977,  1978)  or  wave  2  as  supported  by  Schoeberl  and  Strobe! 
(1980b)  Is  most  important  In  the  warming  process.  The  Cl  warming 
showed  a  greater  wave  2  tendency  with  its  split  vortex  in  contrast  with 
C2,  which  was  more  aligned  with  wave  1.  However,  the  time  sections  for 
both  cases  tend  to  indicate  that  indirect  cells  develop,  westerlies 
decelerate,  and  temperatuers  Increase  most  frequently  when  wave  1  is 
amplifying  and  while  wave  2  Is  decaying.  It  appears  then  that  it  is 
the  push  from  an  Intensifying  wave  which  causes  the  sudden  changes  and 
in  this  context  both  waves  1  and  2  are  Important  since  each  can  serve 


as  an  energy  source  for  the  other  because  of  wave  to  wave  interactions. 
Wave  3  aided  in  this  excahnge  of  energy  serving  as  a  source  of  dissi¬ 
pation  as  proposed  by  Perry  (1967)  and  also  as  an  interloper  linking 
waves  1  and  2. 

In  the  troposphere,  changes  were  not  as  dynamic  as  the  upper 
atmosphere,  probably  due  to  the  anchoring  of  topography  and  trunca¬ 
tion  at  wave  4.  To  some  extent  the  troposphere  did  reflect  the 
stratospheric  warming  event  with  smaller  but  similar  changes  in  zonal 
velocity  and  temperature  at  both  high  and  low  latitudes.  The  main 
differences  between  Cl  and  C2  In  the  troposphere  were  the  Increased 
wave  3  and  4  activity  of  C2,  the  larger  warming  of  polar  regions  that 
occurred  with  Cl,  and  the  easterlies  which  also  appeared  in  the  polar 
regions  of  Cl.  In  both  Cl  and  C2,  there  was  some  tendency  of  tropo¬ 
spheric  geopotentials  to  lead  stratospheric  peaks. 

Finally,  our  results  indicate  the  need  for  the  Inclusion  of  some 
upper  atmospheric  parameterization  to  aid  In  long  range  wintertime 
tropospheric  predictions.  Interactions  between  the  troposphere  and 
stratosphere  resulted  In  enhanced  wave  activity  at  higher  wavenumbers 
which  was  proportional  to  the  strength  of  the  polar  night  jet.  The 
troposphere  also  reflected  the  temperature  and  zonal  wind  changes  that 
occur  during  sudden  warmings  in  the  stratosphere.  The  location  and 
movement  of  troposherlc  jet  maxima  also  tended  to  behave  differently 
and  more  erratically  when  the  upper  atmoshere  contributions  were 
neglected.  The  troposphere  is  not  separated  from  the  stratosphere  in 
nature  and  cannot  be  separated  in  the  modeling  process  if  we  want  to 
accurately  predict  dynamic  meteorological  phenomena. 
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Figure  2.  Latitude-height  sections  of  initial  mean  zonal  wind 


(x  10  m  s'1)  for  (a)  Case  1;  and  (b)  Case  2. 
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Figure  3.  Latitude-height  sections  of  the  initial  zonal  mean 
temperature  delations  (x  10°K)  for  (a)  Case  1;  and 
(b)  Case  2. 
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Same  as  Fig.  5  except  on  (a)  day  22;  (b)  day  24; 
(c)  day  26;  and  {d)  day  28. 


Figure  7.  Latitude-height  sections  of  zonal  mean  temperature 


A 


deviations  (x  10°K)  for  Cl  on  (a)  day  10;  (b)  day 
20;  (c)  day  30;  and  (d)  day  40. 
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Figure  8.  Same  as  Fig.  7  except  on  (a)  day  22;  (b)  day  24; 

(c)  day  26;  and  (d)  day  28. 
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Figure  9.  Cl  latitude-height  sections  of  I *x I  (x  102  m2  s-2)  indi¬ 
cated  by  solid  lines  (except  for  |$x|  =  100  represented 
by  long  dashed  lines)  with  the  phase  angle  (radians) 
indicated  by  short  dashed  lines  for  (a)  day  10;  (b)  day 
20;  (c)  day  30;  and  (d)  day  40. 
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Figure  21 


Cl  latitude-time  sections  at  40.5  km  of  (a)  the 
zonal  mean  meridional  velocity  component  (m  s'1)  and 
(b)  the  zonal  mean  vertical  velocity  (x  10_1  m  s_1). 
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Figure  22 


Cl  latitude-time  sections  of  zonal  mean  temperature 


deviation  (°K)  at  (a)  7.5  km;  (b)  40.5  km;  and 
(c)  67.5  km. 
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Figure  24.  Same  as  Fig.  23  except  for  |$2|. 
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Figure  26.  Same  as  Fig.  23  except  for  I* 
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Figure  30.  Same  as  Fig.  29  except  on  (a)  day  32;  (b)  day  34; 
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(c)  day  36;  and  (d)  day  38. 
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Figure  31.  Latitude-height  sections  of  zonal  mean  temperature 
deviation  (x  10°K)  for  C2  on  (a)  day  10;  (b)  day  20; 
(c)  day  30;  and  (d)  day  40. 
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Figure  32.  Same  as  Fig.  31  except  on  (a)  day  32;  (b)  day  34; 


(c)  day  36;  and  (d)  day  40. 
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Figure  33.  C2  latitude-height  sections  of  I <t> 2 1  (x  102  m2  s-2)  indi¬ 
cated  by  solid  lines  (except  for  1 4^ I  =  100  represented 
by  long  dashed  lines)  with  the  phase  angle  (radians) 
indicated  by  short  dashed  lines  for  (a)  day  10;  (b)  day 
20;  (c)  day  30;  and  (d)  day  40. 
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Figure  34.  Same  as  Fig.  33  except  for  *2. 


Figure  36.  Polar  stereographic  projections  of  height  deviations 
(x  102  m)  for  C2  at  approximately  6.7  km  on  (a)  day 
30  and  (b)  day  40;  and  at  13.5  km  on  (c)  day  34  and 
(d)  day  38. 
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mean  meridional  velocity  component  (m  s_l);  (c)  zonal  mean  vertical  velocit 
(x  10_1  m  s~ 1 ) ;  and  (d)  zonal  mean  temperature  deviation  (°K). 
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Figure  43.  C2  latitude-time  sections  at  40.5  km  of  (a)  the  zonal 
mean  meridional  velocity  component  (m  s"1)  and 
(b)  the  zonal  mean  vertical  velocity  (x  10'1  m  s'1). 
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Figure  46.  Same  as  Fig.  45  except  for  |*2| 
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Figure  47.  Same  as  Fig.  45  except  for  |#  | 
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gure  48.  Same  as  Fig.  45  except  for  | » 
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Figure  55.  C3  latitude-time  sections  at  7.5  km  of  (a)  mean 

zonal  velocity  (m  s_l)  and  (b)  zonal  mean  temperature 
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